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This  dissertation  describes  the  results  of  an  investigation  of 
radiation  transfer  through  aerosols.  Aerosols  have  been  studied  for 
many  years,  and  several  good  reference  works  on  the  physical  properties 
of  aerosols  are  available.  (Green  and  Lane,  1964;  Sedunov,  1974:  Fuchs 
and  Sutugin,  1970;  Junge,  1963;  Voloshchuk  and  Sedunov,  1973;  Friedlander, 
1977.)  Sedunov  presents  a concise,  useful  definition  of  an  aerosol:  a 
solid  and/or  liquid  substrate  in  a gaseous  superstrate.  Basically  an 
aerosol  is  a collection  of  solid  und/or  liquid  particles  suspended 
in  a gas.  1'he  gas  involved  is  usually  an  atmosphere  such  as  that 
of  the  earth,  and  the  particles  are  suspended  due  to  some  act  of 
nature  or  man.  Aerosols  may  be  generated  by  meteorological 
conditions  (condensation  fog),  by  the  wind  (dust  or  near  surface 
haze),  and  by  combustion  (most  smokes). 

Most  liquid  aerosols  are  comprised  of  water,  often  with 
impurities.  Solid  aerosols  are  more  varied,  being  comprised  of  such 
varied  materials  as  minerals  (dust),  elemental  carbon,  or  organic  com- 
pounds (spores).  In  some  situations,  it  is  possible  for  an  aerosol 
particle  to  be  both  solid  and  liquid.  An  example  of  this  would  be  a 
dust  aerosol  that  is  exposed  to  varying  humidity  conditions:  a gust  of 
wind  raises  a cloud  of  dust,  the  particles,  which  are  irregular  in  shape. 
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are  carried  into  a region  of  greater  humidity  and  water  is  adsorbed 
on  the  surface  of  the  particle(s).  This  adsorbtion  has  the  effect 
of  making  the  particle  more  spherical  in  shape  and  increasing  its 
size.  If  the  particle  is  insoluble,  it  may  continue  to  adsorb 
water  until  equilibrium  is  reached.  Sufficient  water  may  be 
adsorbed  to  completely  cover  the  particle.  Then  this  particle 
would  be  layered,  with  an  irregular,  solid  core,  and  coated  with 
liquid  that  is  externally  spherical.  If  the  particle  were  water 
soluble,  the  water  solution  would  eventually  form  a spherical  water 
particle  with  impurities. 

In  general,  liquid  aerosol  particles  are  spherical  until 
they  become  large  enough  for  gravity  and  aerodynamic  drag  to  deform 
them.  Solid  aerosol  particles  are  usually  nonspherical , although 
spores,  some  forms  of  carbon,  and  some  plastics  may  be  spherical. 

Aerosols  are,  of  course,  a collection  of  particles,  and  in 
addition  to  their  shape,  must  be  described  in  terms  of  their  sizes. 
(Aerosol  collections  made  up  of  particles  of  only  one  size  are 
monodisperse,  those  of  particles  of  several  sizes  are  polydisperse. ) 
The  number  of  particles  per  unit  volume  per  unit  size  is  usually 
described  by  a particle  size  distribution  function,  size  usually 
being  expressed  in  terms  of  the  particle  radius,  although  particle 
diameter,  area,  volume  or  mass  may  be  used  (Corn,  1976).  Because  the 
radius  representation  is  the  most  common,  and  is  well  suited  to 
radiation  transfer  calculations,  the  radius  representation  will  be 
used  in  this  dissertation. 
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There  are  three  main  types  of  particle  size  distribution  j 

functions  in  common  use.  These  are  the  Junge  (Junge,  1963), 

n (r)  - C/rn  (1.1) 

where  n(r)  - number  of  particles  per  unit  volume  per  unit  radius, 

r ■ particle  radius, 

C « a constant,  and 
n ■ coefficient, 

the  Deirmendjian  modified  gamma  (DeirmendJ ian,  1964), 

n(r)  - ArVxp  [-Br^j  (1.2) 

where  A,  B,  >,  6,  are  the  parameters,  and 
the  log  normal  (Corn,  1976), 

n(r)  - ^ exp  [-  (ln(r)  - In  (<r*) }2/21n(o)“l  (1.3) 

where  D - coupling  constant, 

*r'  ■ average  particle  radius,  and  , 

o - standard  deviation. 


The  Junge  distribution  is  commonly  used  to  describe  atmospheric  aero- 
sols (haze).  Is  defined  with  an  upper  and  lower  limit  on  particle  radius, 
and  has  a coefficient  n that  is  observed  to  vary  between  the  limits  of 
2.5  and  4.  A Junge  distribution  is  shown  in  Figure  (1.1). 

The  Deirmendjian  distribution  is  commonly  used  to  describe  fog. 
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hate,  and  rain.  An  example  of  the  Delrmendjlan  distribution  is  show 
in  Figure  (1.1).  The  log  normal  distribution  is  also  commonly  used  to 
describe  fogs  and  hates  as  well  as  many  smokes.  This  distribution  enjoys 
much  experimental  popularity  due  to  the  ease  it  affords  data  analysis. 

An  example  of  a log  normal  distribution  is  also  show  in  Figure  (1.1). 

While  nonspherical  aerosols  exist  and  are  commonplace,  they  are 
usually  described  in  terms  of  one  of  these  three  distributions  and  rules 
exist  for  determining  the  "equivalent"  spherical  site  of  the  particles 
(Cadle,  19ns) . No  attempt  is  made  by  Cadle  to  describe  the  shape 
distribution  of  the  particles  as  a function  of  site. 

Aerosols  are  also  described  by  their  refractive  indices.  This 
quantity  is  usually  complex  and  because  of  the  difficulty  of  measuring 
it,  is  assumed  to  be  constant  for  any  one  aerosol  with  respect  to  parti- 
cle site  and  shape.  The  refracttve  index  of  water  is  show  as  an  example 
in  Figure  (1.2)  (Delrmendjlan,  1975). 

The  first  step  in  treating  the  problem  of  radiation  transfer 
through  aerosols  is  to  solve  the  problem  of  the  scattering  of  light  by 

a single  particle.  This  problem  has  been  studied  with  notable  success 
since  the  formulation  of  the  Maxwell  Equations  (Rayleigh,  1871:  Mie, 

1108;  Debye,  19(19).  While  an  analytic  solution  for  the  scattering  and 
absorption  of  a plane  wave  of  light  of  arbitrary  wavelength  by  a sphere 
of  arbitrary  radius  and  refractive  index  (the  Mie  solution)  has  long 
been  know,  the  extent  ot  the  numerical  calculations  has  pre- 
cluded extensive  use  of  the  solution  until  the  recent  advent  ot  the 
digital  computer.  The  Rayleigh  solution  for  particles  small  in 
radius  compared  to  the  wavelength  is  relatively  simple  and  has 
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Figure  (1.2).  Refractive  Index  of  Water 
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therefore  enjoyed  much  wider  use  prior  to  the  computer's  advent.  While 
the  Rayleigh  solution  is  not  applicable  to  most  aerosols  for  wavelengths 
shorter  than  the  submillimeter  (100  pm.),  it  is  applicable  to  gas  molecules 
even  into  the  ultraviolet,  and  much  work  on  the  single  and  multiple  scat- 
tering properties  of  clear  atmospheres  has  been  conducted. 

At  this  time,  no  elegant  solution  to  rival  Mie  theory  has  appeared 
to  describe  the  scattering  properties  of  general  nonspherical  particles 
except  for  small  size  particles  which  may  be  treated  by  Rayleigh-Gans 
theory  (Cans,  1925).  Attempts  have  been  made  to  adjust  particle 
size  distributions  used  with  Mie  theory  to  reproduce  experimental 
data,  but  these  efforts  have  met  with  limited  success  (Cadle,  1965). 

In  recent  years,  several  treatments  of  the  problem  of  single 
scattering  by  irregular  particles  have  been  advanced.  These  treatments 
have  included  extensions  of  the  Mie  solution  into  other  coordinate 
systems  where  Helmholtz's  equation  is  separable.  One  notable  example 
of  this  type  of  effort  is  the  solution  for  orolate  and  oblate 
spnerlods  (Asano  and  Yamamoto,  1975). 

Perturbation  techniques  have  been  applied  by  expanding  the 
electric  and  magnetic  fields,  and  the  particle's  shape  in  perturbation 
series.  Application  of  the  boundary  conditions  results  in  a perturbation 
equation  that  allows  the  corrections  to  the  electric  and  magnetic  fields 
to  be  calculated  in  a bootstrap  manner  (Yeh,  1964).  Unfortunately, 
this  method  suffers  from  the  usual  limitation  of  perturbation  theory  — 
the  perturbations  must  be  small. 

A technique  that  has  enjoyed  some  use  in  describing  conductors 
is  the  use  of  the  last-squares  regression  technique  to  solve  for  the 
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scattered  electric  and  magnetic  fields  (Davies,  19731.  This  technique, 
however,  does  not  exactly  satisfy  the  boundary  conditions  of 
e lec  t rodynamics . 

One  recent  technique  of  note  has  been  the  use  of  an  integral 
equation  formalism  (Eyres  and  Nelson,  1976;  Uzonoglu  and  Holt,  1977; 
Waterman,  1971).  While  Integral  equation  formalisms  have  been  used 
before,  these  latest  treatments  are  notable  because  the  boundary  con- 
ditions are  properly  accounted  for,  and  the  numerical  solution  Is 
obtained,  not  by  iteration,  but  by  reducing  the  Integral  equations  to 
linear  equations  and  obtaining  the  expansion  coefficients. 

An  empirical  treatment  recently  advanced  is  a modification  of 
the  Mie  solution  based  on  the  correlation  of  resonances  in  the  ana- 
lytical solution  with  glories  in  the  experimental  data  (Chylek,  Grams, 
and  Pinnick,  1976),  Although  calculations  performed  using  this  modi- 
fication have  displayed  better  agreement  with  experimental  data  than  Mie 
calculations,  the  validity  of  the  modification  has  been  criticized 
(Kerker,  1977).  While  this  treatment  offers  a possibility  of  widespread 
usefulness  because  of  its  computational  simplicity,  the  modification's 
basis  has  not  been  substantiated  by  more  exact  calculations. 

The  second  step  in  treating  the  radiation  transfer  problem  Is 
to  extend  the  calculations  of  the  single  scattering  of  light  by  a 
single  particle  to  the  single  scattering  of  light  by  a collection  of 
particles.  In  most  cases,  notably  for  naturally  occuring  aerosols  and 
the  visible  wavelengths  of  light,  aerosol  particles  are  sufficiently 
large  that  the  Mie  solution  Is  required.  Fortunately,  the  aerosols  are 
also  sufficiently  dilute  so  that  during  any  reasonably  long  period  oi 
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time,  the  particles  are  in  the  far  fields  of  all  the  other  particles. 
DeirmendJ ian  has  shown  that  under  these  circumstances,  the  particles 
are  independent  and  the  scattered  intensities  rather  than  the  scattered 
electric  and  magnetic  fields  maj  be  added  (Diermendj ian,  1964).  As  a 
result,  the  scattering  properties  of  most  aerosols  may  be  described  in 
terms  of  an  average  particle.  It  should  be  noted,  however,  that  work 
has  been  done  to  include  the  interparticle  correlations  (Foldy,  1945). 

As  the  submillimeter  and  millimeter  regions  of  the  spectrum  become  more 
accessable  and  utilized,  the  importance  of  this  problem  should  become 
more  apparent. 

The  third  step  in  treating  the  problem  of  radiation  transfer 
through  aerosols  is  the  calculation  of  the  multiple  scattering  effects 
that  may  occur  in  an  aerosol  medium.  To  perform  this  calculation,  the 
theory  of  radiative  transfer  (RT)  developed  by  Chadrasekhar  may  be 
used  (Chandrasekhar,  1960).  This  theory  has  enjoyed  extensive  use  in 
describing  the  propagation  of  light  through  stellar  and  planetary 
atmospheres.  These  two  problems  are  characterized  by  three  notable 
conditions:  first,  the  illumination  incident  onto  and/or  from  the  atmos- 
phere is  uniform;  second,  the  scattering  medium  varies  in  only  one  direc- 
tion at  most;  and  third,  the  scattering  phase  function  of  a gas  (the 
Rayleigh  solution)  has  a simple  analytic  form.  These  conditions  allowed 
the  general  form  of  the  RT  equation  to  be  reduced  to  a one-dimensional 
form  that  could  in  a few  cases  be  solved  analytically  and,  for  the 
Rayleigh  phase  function,  could  be  solved  numerically  without  the  computer. 
Indeed,  this  numerical  solution  represented  the  essential  state-of-the- 
art  prior  to  the  computer's  advent  about  1960. 
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The  digital  computer  has  permitted  the  RT  equation  in  its  one- 
dimensional form  to  be  solved  for  other,  more  complicated  phase  func- 
tions so  that  the  problems  of  light  transmission  through  a hazy  atmosphere 
or  an  atmosphere  with  uniform  cloud  cover  could  be  solved. 

While  this  one-dimensional  form  of  the  RT  equation  was  adequate 
for  describing  the  brightness  of  a star,  or  the  amount  of  light  reaching 
the  surface  of  a planet,  advances  in  detector /source  technology  have 
uncovered  problems  that  the  existing  RT  techniques  cannot  treat.  The 
availability  of  both  visible  and  infrared  imaging  devices  and  of  the 
laser  have  introduced  a new  class  of  problems  involving  nonuniform 
incident  Illumination.  The  determination  of  laser  beam  broadening  or 
imaging  contrast  degradation  due  to  multiple  scattering  is  beyond  the 
scope  of  the  uniform  illumination  formulation  of  RT. 

To  date,  however,  only  limited  attempts  have  been  made  to  solve 
the  general  RT  equation.  Notable  among  these  is  a calculation  using  the 
four-dimensional  RT  equation  simplified  for  small  angle  scattering  only 
(Weinman  and  Shipley,  1972).  This  calculation,  which  was  performed  to 
predict  laser  pulse  stretching  in  clouds,  is  valid  only  when  the  aerosol 
particles  are  much  larger  than  the  wavelength.  (The  phase  function  is 
sharply  peaked.)  This  approach  is  not  valid  for  scattering  by  atmos- 
pheric gases  because  the  Rayleigh  phase  function  is  not  sharply  peaked. 

Another  notable  calculation  was  performed  to  support  the  optical 
measurement  of  the  particle  size  distribution  in  a chamber  (Deepak  and 
Green,  1975).  Although  this  calculation  was  valid  for  general  scatter- 
ing angle,  it  was  carried  out  only  to  third  order  in  scattering.  While 
third  order  scattering  is  adequate  for  the  optical  depths  observed  in 
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the  experiment,  the  technique  of  Deepak  and  Green  is  not  readily  applic- 
able to  problems  that  require  that  higher  orders  of  scattering  and/or 
more  general  boundary  conditions  be  considered. 

As  a result  of  these  considerations,  this  dissertation 
investigation  was  begun  to  determine  the  effects  of  aerosol  and 
atmospheric  gas  multiple  scattering  on  imaging  and  laser  beam  trans- 
mission. These  transmission  calculations  involve  nonuniform  incident 
illuminations  and  aerosol  media  that  had  not,  to  date,  been  generally 
treated,  either  analytically  or  numerically.  This  phase  of  investi- 
gation has  resulted  in  a two-dimensional  RT  code  for  the  treatment 
of  imaging  and  laser  beam  transmission.  This  code  constitutes  the 
contribution  of  this  investigation  to  the  third  step  of  the  radiation 
transfer  problem. 

As  calculations  with  this  code  proceeded  from  liquid  aerosols  to 
solid  aerosols,  the  problem  of  single  scattering  by  irregular  particles 
surfaced.  As  has  been  noted  earlier,  several  techniques  for  treating 
the  problem  were  known,  but  none  provided  the  necessary  information  for 
solution  of  the  RT  equation.  This  information  includes  the  extinction 
coefficient,  the  albedo,  and  the  phase  functions  for  irregularly  shaped 
dielectric  aerosols.  Considerable  effort  has  been  made  to  extend  the 
Mie  solution  to  irregular  particles  to  provide  this  information.  This 
extension  was  implemented  into  a code  for  solving  the  single  scattering 
problem  for  cylindrically  symmetric  irregular  dielectric  particles.  This 
restriction  is  imposed,  not  by  limitations  of  the  theory,  but  by  limi- 
tations of  the  computer.  The  development,  implementation  and  use  of 
this  code  is  described  in  this  dissertation. 
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This  code  constitutes  the  contribution  of  this  investigation  to 


the  first  part  of  the  problem  of  radiation  transfer  through  aerosols. 


Thus,  of  three  parts  to  the  problem,  two  have  been  addressed  in  this 


investigation:  the  single  scattering  by  an  Irregular  dielectric  particle. 


and  the  multiple  scattering  effects  for  nonuniform  boundary  conditions. 


In  Chapter  II  of  this  dissertation,  the  Hertz  vector  formalism 


is  developed  from  the  Maxwell  Equations.  This  formalism  has  the  ad- 


vantage of  reducing  the  vector  differential  equations  (which  are  derived 


directly  from  the  Maxwell  Equations)  for  describing  the  single  scatter- 


ing to  scalar  differential  equations.  This  Hertz  vector  formalism  is 


then  used  in  Chapter  III  to  derive  the  Mie  solution.  The  extension  of 


the  Mie  solution  to  aerosols  is  reviewed  in  Chapter  IV. 


Several  of  the  previous  formalisms  for  treating  scattering  by 


irregular  particles  are  reviewed  and  discussed  in  Chapter  V.  Chapter  VI 


then  presents  the  theory  for  single  scattering  by  irregular  particles 


used  in  this  investigation.  This  is  an  extension  of  the  differential 


equation,  Hertz  vector  formalism  of  Mie  to  "slightly"  irregular  particles 


with  cylindrical  symmetry.  Unlike  the  Mie  solution  for  the  spherical 


particle,  which  has  an  analytic  solution  that  can  be  evaluated,  this 


formalism  requires  a numerical  solution  because  the  radius  of  the 


particle  is  itself  a function  of  angle.  This  solution  may  be  obtained 


by  the  computer,  however,  and  a code  to  perform  this  calculation  is 


described  in  Appendix  (I). 


The  calculations  performed  using  this  code  for  spheres  and 


spheroids  are  compared  with  exact  calculations  in  Chapter  VI  l.  Addi- 


tionally, sample  calculations  for  a variety  of  irregular  particles  are 


Ift 


I 


Y I 


y 


presented.  In  Chapter  VIII,  the  validity  of  the  Mie  modification  of 
Chylek,  Grams  and  Pinnick  is  addressed  (Chylek,  Grams,  and  Pinnick, 
1976).  This  modification  is  first  examined  by  calculation  of  the 
electric  and  magnetic  field  expansion  coefficients.  Following  this, 
polydisperse  aerosol  phase  functions  are  calculated  for  several  types 
of  irregular  particles  after  incorporating  the  techniques  of  Chapter  IV 
into  the  code.  These  phase  functions  are  then  combined  using  linear 
regression  techniques  to  fit  experimental  data. 

Chapter  IX  begins  the  discussion  of  the  multiple  scattering  part 
of  this  investigation.  In  this  chapter,  the  general  RT  equation  is  pre- 
sented and  its  reduction  to  one-dimensional  form  for  the  treatment  of 
plane-parallel  atmospheres  is  outlined.  An  approximate  solution  of  the 
two-dimensional  unpolarized  RT  equation  suitable  for  imaging  and  steady- 
state  laser  beam  transmission  calculations  is  presented  in  Chapter  X, 
together  with  an  algorithm  for  its  numerical  solution.  This  algorithm 
is  implemented  in  a code  described  in  Appendix  (IV).  Calculations  per- 
formed for  transmission  through  a Deirmendjian  C.3  fog  using  this  code 
are  presented  in  Chapter  XI  and  discussed.  Chapter  XII  is  a summary  of 
the  results  of  this  dissertation  investigation  with  emphasis  on  the 
various  limitations  of  the  calculations.  A few  problems  to  be  treated 
as  future  efforts  are  indicated. 
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CHAPTER  II 


ELECTROMAGNETIC  WAVES  AND  HERTZ  VECTORS 


derived  from  Maxwell's  Equations  (Jackson,  1962), t 


the  electric  displacement 


the  magnetic  flux  intensity 


the  magnetic  field 


the  speed  of  light 


The  constitutive  equations  are 


t Gaussian  units  are  used  throughout  this  dissertation 


H . (2.13) 

This  set  of  Maxwell  Equations,  Eqs.  (2.8),  (2.9),  (2.12),  and  (2.13)  may 
be  used  to  calculate  the  light  scattering  properties  of  individual 
particles . 

In  particular.  Equations  (2.8)  and  (2.9)  admit  the  solution  of 
the  Maxwell  Equations  by  Hertz  vector.  Because  the  divergence  of  the 
curl  of  a vector  is  zero  (V*VxA  - 0),  the  magnetic  and  electric  Hertz 
vectors  may  be  defined  as 

| - aV*|l,  (2.14) 

and 

IJ  - bVx£  , (2.15) 

where  a and  b are  coefficients  to  be  determined.  It  may  be  noted  that 
]\  and  ^ are  two  Independent  solutions  of  the  second  order  vector 
partial  differential  equations  derivable  from  the  Maxwell  Equations. 
Because  ^ and  ^1  must  be  transverse  in  the  far  field,  J1  and  £ must  be 
radial  vector  functions  and  any  linear  differential  equations  in  J1 
and  £ must  be  reducable  to  scalar  equations. 

The  magnetic  field  due  to  J^I  may  be  found  by  substituting 
Eq.  (2.14)  into  Eqs.  (2.12)  and  (2.13)  to  yield 

Vxj(  - - tkam2V*jl  , (2.16) 

and 

• (2*17) 
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so  Chat 


- • V . ”r.  * 





•^™r- 


^ ■ -tkam2]!  + V<t>  , (2.18) 

where  Hs  a scalar  potential  that  Is  consistent  with  Eq.  (2.16)  since 
the  curl  of  the  gradient  of  a scalar  is  zero.  This  potential  is 
included  to  insure  that  jl  will  satisfy  the  Lorentz  condition  (Jackson, 


1962;  Tyras,  1969), 

* - . (2.19) 

Equation  (2.17)  may  be  rewritten  using  the  vector  identity, 

z*i*b  - w$>  * (2-20) 

to  become 

- a^n  “ *klJ  » (2<21) 

which  may  be  combined  with  Eq.  (2.18)  to  yield 

jW(V-JT)  - aV2Jl  - ak2m2JT  + ^.kV*.  (2.22) 


Equation  (2.22)  may  be  combined  with  Eq.  (2.19),  the  Lorentz  condition, 
to  yield 

aV(V*n)  - aV2JJ  - ak2m2JJ  + ikV(V«£)  . (2.23) 

Since  a is  an  arbitrary  constant,  it  may  be  identified  as  ^tk,  so  that 
Eq.  (2.23)  reduces  to 

V2]l  + kJm2]l  - 0 , (2.24) 

Helmholtz's  wave  equation.  Equations  (2.14)  and  (2.18)  may  be  rewritten 
as 


I, 
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I - <kW  • 

and 

« “ ^ + WJP 


(2.25) 


(2.26) 


The  electric  Hertz  vector,  £,  may  be  derived  in  a similar 
manner,  starting  with  the  substitution  of  Eq.  (2.15)  into  Eqs.  (2.12) 
and  (2.13).  This  substitution  yields 

bV*Vx£  « - <km^|  (2.27) 


and 

Vx£  - ikbV*£  . (2.28) 

Equation  (2.28)  may  be  exploited  to  admit 

£ - tkb£  + VY  , (2.29) 

where  Y Is  a scalar  potential  that  satisfies  the  Lorentz  condition 

Y ■ . (2.30) 

Equation  (2.27)  may  be  rewritten  using  Eq.  (2.20)  as 

bV(V*£)  - bV2*  m bk2m2£  - -ikm^Y  . (2.31) 

, 2 

Equation  (2.27)  allows  b to  be  identified  as  -tkm  , and  Eq.  (2.31) 
becomes 

V2£  + k2m2£  - 0 , (2.32) 

the  Helmholtz's  wave  equation  for  the  Magnetic  Hertz  vector.  Further, 
Equations  (2.15)  and  (2.29)  may  be  rewritten  as 

- -<km2V*J[  , (2.33) 
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which  have  general  solutions  of  the  forms  (Born  and  Wolf,  1975) 


(cost)) 


a»i> » (krar ) + b,,n,>(kmr) 


(kmr)  + bap^(kmr)  PT(cos0) 


depending  on  the  boundary  conditions  to  be  applied.  The  angular 
function  is  the  associated  Legend  polynomial  of  order  i. 


are  related  to  the  usual 


the  cylinder  functions  tK,  n 


spherical  Bessel  functions  (Abramowitz  and  Stegun,  1964)  by 


Equations  (2.37)  and  (2.38)  may  be  explicitly  written  as 
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CHAPTER  III 


SPHERICAL  p ARTICLES :MIE  THEORY 


The  Mte  problem  Is  that  of  describing  the  scattered  electric  and 
magnetic  fields  due  to  the  interaction  of  a plane  wave  of  light  with  a 
spherical  object.  The  boundary  conditions  for  this  problem  require  that 
the  components  of  the  fields  tangent  to  the  sphere  at  its  surface  be 
continuous.  Because  the  surface  normal  for  a sphere  is  purely  radial, 
the  boundary  conditions  require  that  E(,,  E^,  H(1,  and  be  continuous 
at  the  sphere's  surface.  Examination  of  Eqs.  (2.49)  - (2.55)  reveals 
that  if  n and  o are  continuous  and  once  radially  differentiable  contin- 
uous at  the  surface,  the  boundary  conditions  will  be  satisfied  since 
the  angular  dependence  of  Equations  (2.43)  and  (2.44)  are  Identical. 

The  plane  wave  is  taken  to  be  Incident  along  the  z axis  with 


components 


r *kz 
L “ e , 
x 


(3.1) 


Hy  " 


(3.2) 


which  give  rise  to  scalar  potentials  (in  spherical  coordinates)  of 


i -2  v./-l  2(  + 1 
0 “ k Y 1(1  + ! 


(kr)  pj.  (cos  0)  cos  ( $) , (3.3) 


n1  - k"2  £ t*"1  j-^VlT  *t  (kr)  Pf  (cOS  0)  8ln  U)*  (3,4) 


1 


A - . - _ J ' 


where  the  superscript  "i"  denotes  the  incoming  wave  (solution)  and  m ■*  1 


(free  space)  outside  the  sphere.  The  assymptotic  behavior  of  the 


cylinder  functions  in  the  limits  r ■+■  0,  and  r -*■  <»  allow  the  internal 


("w")  and  scattered  ("s")  potentials  to  be  written  as 


ow  k m (kmr)  (cos  6)  cos  (<p) , 


(3.5) 


k-2  m 1 (kmr)  (cos  0)  sin  (<J>), 


(3.6) 


os  - k-2  (kr)  P^  (cos  0)  cos  (4>)  , 


(3.7) 


s , -2 


k_2  (kr)  P^  (cos  0)  sin  ($) 


(3.8) 


The  difference  of  a factor  m * between  Eqs.  (3.5)  and  (3.6)  arises  from 


the  fact  that  if  £ = e^  Ae^™2  and  ^ Be^0112,  then  Eqs.  (2.12)  and 


(2.13)  require  that  B ■ mA.  The  boundary  conditions,  that  the  tan- 


gential components  of  the  £ and  fields  be  continuous,  may  be  determined 
by  examination  of  Eqs.  (2.49)  - (2.54).  Since  v and  o are  already 


continuous  in  angles,  the  tangential  continuity  is  satisfied  by:  Eg, 


and  ; E , tt  and  |~;  H0,  and  m2o;  and  H^,  and  m o.  Thus,  the 


boundary  conditions  are  specifically. 


— (a^  + os)  | * t—  ow| 

r v 'r  - a 9r  ‘r  - a. 


9 _wi 


(3.9) 


3 / 1 S\  I wi 

9r  n ^ ^r  - a ” 9r  U • a. 


9 _wi 


(3.10) 


(o1  + os) | r 


2 w | 

m o 

'r  ■ a. 


(3.11) 


m 


where  a Is  the  radius  of  the  sphere 


The  orthogonality  properties  of  sin(<J>),  cos(<j>),  and  the 


(cos *P)  may  be  exploited  to  reduce  Eqs.  (3.3)  - (3.12)  to 


where:  tp' (y) 


Equations  (3.13)  - (3.16)  may  be  combined  to  solve  for  the  expansion 


coefficients  of  the  scattered  wave  scalar  potentials 


(kma)p 


By  combining  Eqs.  (2.50),  (2.51),  (2.53),  (2.54),  (3.7),  (3.8), 
(3.17),  and  (3.18),  the  scattered  electric  and  magnetic  field  and 
components  may  be  written  as 


Ee  - (kr) 


-1 


, Pj(cose)  3P^1(cos  0) 

P£  (kr)  sln(Q)  + c£p£  (kr)  53- 


X cos  ($)  , (3.19) 


E - 


•(kr)'1^ 


. 3P»(cos0) 

*** lPl  *kr*  36 


, P„  (cos  0) 

cA  (kr)^n?) 


X sin  (<#>) , (3.20) 


H0  " )~ll 


, 3Pt(cos0) 

V?  <*>  sf 


. P (cos  0) 

- ictpi  (kr)  Ye") 


X sin  (4>)  , (3.21) 


<kr)_1£ 


, (CO80) 

dp1  Otr)  — 


l l 


sin  (9) 


. 3P»  (cos  0) 

+ ictP(  (kr)  Tq 


Xcos($).  (3.22) 


The  radial  components  may  be  calculated  in  the  same  manner,  but  are 

-2  -1 

not  included  because  they  decrease  as  r rather  than  as  r and  will 
be  vanishingly  small  in  the  far  field  when  compared  to  Eqs.  (3.19)  - 


(3.22). 


The 


asymptotic  properties  of  the  p»*  's  are 
p£l (x)  - (- i)t+1  e**,  x - - 


(3.23) 


x 00 


(3.24) 


O’ 
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which  may  be  used  to  write  the  far  field  components  as 


; 1 1 1 i ■■  i i . 


Ee  " 


^kr  p. , 

e r / • v c+1 

It  (“O 


kr 


l 


f ?i^co 

rl  sin(0) 


cos  0)  3P»  (cos  0) 

- + c 


'l  30 


cos  (4.), 


-ckr 
e 

kr 


(-*) 


3P^  (cos  0)  P^(cos  0)  ) 

do  — + c — 


*t  30 


i sin(Q) 


ckr 


'6  kr 


<-<>t+l 


SP^1  (cos  0)  P^1  (cos  0)) 


30 


+ c 


t sin  (0) 


’4>  kr 


4-kr  o. . 


P^(cos  0)  3P^(co8  0) 

dl  sin  (0)  + Cl  30 


It  may  be  noted  that 

E - H 

0 <P 


and 


E*  “ -V 


(3.25) 


sin  (40 , 


(3.26) 


sin  (40 , 


(3.27) 


sin  (40. 


(3.28) 


(3.29) 


(3.30) 


Equations  (3.25)  and  (3.26)  may  conveniently  be  rewritten  (Van  de  Hulst, 
1957)  as 


-ckr 


Ee  ■ • ^ S1(0)  cos(4>). 


(3.31) 


29 
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if 


■ 


and 


kr 

^ " - fe—  S2(0)  sin(<^) 


(3.32) 


where 


S1(6) 


V / ■ \ 

I 


' P„1(CO80)  1 


+ C, 


3P^  (COS0) 


Al  sin(0)  l 30 


(3.33) 


and 


S2(0) 


(-<-) 


l+l 


P (cos0)  3P£1(cos0) 

c o . _ / q\  d 


sin(0)  u£  30 


(3.34) 


A new  set  of  angular  functions  are  commonly  defined  (Van  de  Hulst,  1957; 
Deirmendjlan,  1964), 


P » (cos0) 

Ve)  - 


(3.35) 


and 


" 30 


3P^  (cos0 ) 


(3.36) 


It  is  also  convenient  to  define  new  coefficients  a ^ and  3^, 


ac 


(-<-) 


t+ 1 


c 


til  + 1) 
l 21  + 1 


(3.37) 
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ipo(kma)iK(ka)  - miK(kma)iJiUka) 


(kma)pp(ka)  - miK  (kma)p 


m^Ukma)^p(ka)  - <K(kma)i|>p(ka) 


mi^p(kina)pp(ka)  - ijip(kma'p 


Equations  (3.33)  and  (3.34)  may  now  be  rewritten  as 


It  Is  now  possible  to  develop  the  cross  sections  and  phase 


functions  by  Introducing  the  Stokes  vector  formalism  (Van  de  Hulst 


1957).  The  components  of  the  Stokes  vector  are  defined  as 


where  the  bar  indicates  an  average  over  polarization  angle,  and  E»  and 


E are  field  components.  For  the  spherical  (Mie)  problem  treated  here, 

a 

- -E^,  and  E^  - E0,  so  that 


* 2 „ 
2 S1S1cos  (^) 

h " 7T5 

k r 


(3.47) 


* 2 . 
2 SjSjCos  ($) 

E*  “ , 2 2 
k r 


(3.48) 


EtE-  * JT 

k r 


(S^S*  + S*S2)sin($)cos($) 


(3.49) 


<(S1S2  " sJS2)sin(^)cos(^) 

'*■ ' zn 

k r 


(3.50) 


By  averaging  over  4%  the  differential  cross  sections  may  be  derived  as 

* 


S1S1 

o (a,  X,  m,  0)  - — r — , 
1 k 


(3.51) 


S2S2 

o.(a,  X,  m,  0)  - —z — 
k 


(3.52) 


o3(a,  X,  m,  0) 


(SjS*  + s*s2) 


(3.53) 


X(s1s2  - sjs2) 

2k2 


o^U,  X,  m,  0) 


(3.54) 


""T-T ■»'  — 


' 

I 1 


I | 


where  the  cross  sections  o^  have  been  written  with  all  explicit 
dependence  on  particle  radius,  wavelength,  refractive  index,  and 
scattering  angle.  The  cross  sections  are  related  to  the  Stokes  vector 
components  by  the  matrix  equation 


(3.55) 


l°l 

0 

0 

In o\ 

■ 

f 

0 

°2 

0 

0 

■to 

u 

0 

0 

°3 

°4 

U 

o 

'1 

1° 

0 

-°4 

N 

where  the  subscript  "o"  indicates  the  pre-scattering  value. 
The  total  scattering  cross  section  is  given  by 

osc(a,  X,  m)  - / (ci1  + o2) 

4n 


(3.56) 


where  the  Integral  is  over  solid  angle.  Equation  (3.56)  may  be 
integrated  exactly  for  the  Mie  problem  because  of  the  orthogonality 
properties  of  the  Legendre  polynomials,  yielding 


sc 


(a,  X,  m)  - l (2f  + D^o^,  + 


(3.57) 


The  extinction  cross  section  may  be  evaluated  using  the  optical  theorem 
(Van  de  Hulst,  1957)  as 


oex(a,  X,  m)  - ~ Re  J S^O)  + S2(0)  J 


(3.58) 


The  two  angular  functions  n^(6)  and  1^(8)  may  be  evaluated  at  6 ■ 0 as 
nf(0)  - tf(0)  - +--1*  , 
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(3.59) 


reducing  Eqs.  (3.41)  and  (3.42)  to 

S1(0)  - S2(0)  (2*  + l)Ja^  + 6^]. 

Equation  (3.58)  then  becomes 

oex(a,  X,  m)  - ^ l (2 t + 1)  (Re(o^)  + Re(6^)}  . 
k l 


(3.60) 


(3.61) 


These  equations  may  be  used  to  calculate  the  extinction  properties  of 
the  entire  medium  of  aerosols  as  will  be  discussed  in  the  next  chapter. 


CHAPTER  IV 

SINGLE  SCATTERING  FOR  MANY  SPHERICAL  PARTICLES 
Chapter  III  discussed  the  scattering  of  light  of  arbitrary 
wavelength  by  a single  spherical  particle  of  arbitrary  radius  a and 
refractive  index  m.  In  nature,  aerosol  particles  are  generally 
distributed  in  size  (radius)  although  the  particles  are  usually  either 
spherical  as  in  the  case  of  fog,  or  are  treated  as  spherical  as  in 
the  case  of  atmospheric  haze  or  dust.  Some  of  the  aspects  of  light 
scattering  by  nonspherical  particles  will  be  discussed  in  subsequent 


chapters. 

In  this  chapter,  the  aerosol  particles  will  not  only  be  assumed 
to  be  spherical,  but  also  each  particle  will  be  in  the  far  field  of 
all  other  particles.  If  d is  the  average  distance  between  particles, 
then  it  is  required  that  X<<d,  a<<d.  Since  the  aerosol  particles  are 
moving,  it  is  assumed  that  the  average  time  between  collisions  is 
sufficiently  large  that  the  number  of  particles  for  which  X£d  at  any 
Instant  is  small  with  respect  to  the  total  number  of  particles.  These 
restrictions  allow  the  particles  to  be  treated  as  independent  and  the 
resulting  intensity  quantities,  rather  than  the  field  amplitudes,  are 
summed  over  all  particles. 


It  will  be  further  assumed  that  there  exists  a usable  analytic 

or  numeric  (piecewise  analytic)  function  N(^,  a,  t)  such  that 

3 

N(^,  a,  t)  d r da  dt  is  the  number  of  particles  of  radius  between  a 

3 

and  a+da  within  the  volume  d ^ during  the  time  interval  between  t and 
t+dt.  While  this  function  may  in  principle  be  calculated,  this 
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calculation  has  never  been  performed  exactly,  although  several 
approximate  calculations  have  been  performed  (Voloshchuk  and  Sedunov, 
1973;  Sedunov,  1974). 

While  the  refractive  index  for  any  one  particle  is  assumed 
constant,  the  refractive  index  may  vary  in  position,  time,  or  particle 
size,  corresponding  to  such  cases  as  an  evolving  acid  fog  or  a mixture 
of  water  haze  and  dust  haze.  This  variation  is  represented  by  a 
function  m(^,  a,  t)  representing  the  refractive  index  of  particles 
of  size  a at  position  ^ at  time  t.  Because  of  the  difficulties  in 
measuring  refractive  index,  most  calculations  are  performed  with  a 
constant  refractive  index.  Only  in  some  special  cases  when  it  is 
known  that  two  types  of  aerosol  are  mixed,  such  as  from  a gun  blast 
on  land,  is  the  refractive  index  varied  in  practice. 

Further,  most  experimental  data  on  the  particle  distribution 
function  is  limited  to  one '.(Qr  .at  .best  .a  few) ' measurements  ;o£  the- 
size  distribution  at  a point  averaged  over  some  measurement  time  and 
the  total  concentration  of  aerosols  at. several  points.  As  a result, 
the  particle  distribution  function  N(^,  a,  t)  is  commonly  factored 

N(£,  a,  t)  = f(£,  t)n(a) , (4.1) 

where  f(^,  t)  is  the  concentration  of  particles  and  n(a)  is  called 
the  particle  size  distribution.  While  Eq.  (4.1)  is  known  to  be  true 
for  only  limited,  special  cases,  it  is  commonly  used  as  an  approxi- 
mation for  N(^,  a,  t). 
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Given  N(^,  a,  t)  and  m(^r,  a,  t),  the  extinction  coefficient 
per  unit  length  is  given  by 


a(£,  t)  = /a  (a,  A*  m^»  a»  a»  t)  da* 


(4.2) 


and  the  scattering  coefficient  per  unit  length  is 


tJ(^,  t)  = /o  (a,  X,  m(£,  a,  t))N(jr,  a,  t)  da. 


(4.3) 


where  0 (a,  X,  m)  and  a (a,  X,  m)  were  defined  in  Eqs.  (3.57) 

sc  ex 

and  (3.61).  The  scattering  phase  functions  are  defined  by 


Pi  (£»  t.  9)  = l°.(a>  X,  m(j)r,  a,  t),  0)N(r,  a,  t)  da,  (4.4) 


where  the  Oj(a,  X,  m,  0)  were  defined  in  Eqs.  (3.51)  - (3.54). 
unpolarized  light  phase  function  is  defined  by 


oo 

P(£,  t,  e)  = e(£,  t)-1  /°sc(a.  X,  n>(^»  a,  t)),  0) 


(4.5) 


X N(]£,  a,  t)  da. 


It  may  be  noted  that  P (^ , t,  0)  has  unit  integral  over  solid  angle. 

It  is  convenient  to  define  the  albedo  of  single  scattering  as 


u>(jr,  t)  = 6(£,  t)/a(£,  t), 


(4.6) 


h •£*'  i 
i •*  . 
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As  noted  previously,  Eq.  (4.1),  the  particle  distribution  function 

Is  often  factored  into  two  parts,  f(r,  t)  and  n(a).  The  concentration 

% 

may  frequently  be  represented  by  either  a uniform  cloud,  either  of 
infinite  extent  (fog)  or  of  finite  extent  (patchy  fog),  or  by  some 
concentration  function  such  as  a gaussian  plume  (smoke).  The  particle 
size  distribution  Is  most  commonly  represented  by  one  of  three 
functions;  the  Junge  distribution,  the  log  normal  distribution,  or 
the  Deirmendjiau  modified  gamma  distribution.  The  Junge  distribution 
has  the  form 


n(a)  - C/a  , 


(4.7) 


where  C is  an' adjustable  coefficient  to  reconcile  the  concentration, 
and  n is  usually  of  the  order  of  3 for  most  atmospheric  hazes  (Junge, 
1963).  This  type  of  distribution  is  most  commonly  used  to  describe 
atmospheric  hazes.  The  log  normal  distribution  has  the  form 


n(a)  - ^-exp[-(  (ln(a)  - ln(<r>))/ln(o)  )2/2] , 


(4.8) 


where  D is  again  a coefficient  to  reconcile  the  concentration,  <rv 
Is  the  average  particle  radius,  and  o i«  the  standard  deviation 
(Green  and  Lane,  1964).  This  distribution  is  commonly  used  to  describe 
fogs  and  smokes.  The  Deirmendjian  modified  gamma  distribution  has 
the  form 


Y <5 

n(a)  • Aa  exp (-Ba  ) , 


(4.9) 
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where  A Is  again  a coefficient  to  reconcile  the  concentration  and  B, 
y,  and  6 are  parameters  to  characterise  the  aerosol.  This  distribution 
has  been  used  to  describe  fog,  haze,  snow,  rain,  clouds,  dust,  ami 
smoke  (DeirmendJ tan,  1964).  In  most  cases,  the  parameters  n»  <r>,  o, 

B,  y , and  6 are  measured  experimentally.  While  the  log  normal  and 
DeirmendJ ian  distributions  may  be  reconciled  with  the  overall  concen- 
tration by  a direct  integration,  this  reconciliation  for  the  Junge 
distribution  requires  the  specification  of  the  minimum  radius  to 
avoid  the  pole.  (This  lower  size  limit  also  applied  to  Eqs.  (4.2)  - 


(4.5), 


When  N(£,  a,  t)  is  separable,  and  m is  a constant,  Eqs.  (4.2) 


and  (4.3)  reduce  to 


w 

«(£,  t)  - f(£,  t)joex(a,  X,  m)n(a)da, 


(4.10) 


0($,  t)  - f(£,  t)jogc(a,  X,  m)n(a)da. 


(4.11) 


and  u)(^,  t)  and  P(^,  t,  0)  become  constant  with  respect  to  ^ and  t. 

In  this  case,  the  notation  u)  and  P(0)  is  commonly  used  for  the  albedo 
and  unpolarlzed  phase  function. 
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When  N(^,  a,  t)  Is  not  separable,  an  effective  albedo 
that  Is  the  maximum  value  of  u)(%,  t)  is  commonly  defined.  Since 
u)(^,  t).<l,  u>*<1.  A modified  phase  function 

<*>(£,  t)P(^,  t,  9) 

Q(^,  t,  6)  = — — * (4.12) 

0) 

is  defined  that  has  integral  value  -1  for  all  t.  The  effective 
albedo  and  modified  unpolarized  phase  function  are  commonly  used  in 
Radiative  Transfer  to  permit  an  order-of-scattering  analysis  to  be 
performed  even  when  the  albedo  is  not  constant. 
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CHAPTER  V 

PREVIOUS  STUDIES  OF  NONSPHERICAL  PARTICLE  SCATTERING 
While  the  solution  for  the  scattering  of  a plane  wave  of  light 
of  arbitrary  wavelength  by  a sphere  of  arbitrary  radius  and  refractive 
index  has  long  been  known,  (Mie,  1908;  Debye,  1909)  no  such  elegant 
formalism  has  existed  for  nonspherical  particles  of  arbitrary  size 
and  shape.  There  are  only  a few  geometric  shapes  where  the  Mie 
problem  can  be  solved  exactly  by  using  the  standard  separation  of 
variables  method  (Jackson,  1962).  An  excellent  example  of  this  type 
of  effort  is  the  recent  work  on  prolate  and  oblate  spheroids  (Asano  and 
Yamamoto,  1975).  It  must  be  noted,  however,  that  even  when  the 
analytical  solution  is  available,  the  results  are  extremely  complicated. 

Particles  that  are  small  may  be  treated  with  Rayleigh-Gans 
theory  (Rayleigh,  1881;  Gans,  1925).  This  theory  requires  that 

|m  - 1|  <<1,  (5.1) 

the  refractive  index  of  the  particle  be  only  slightly  different  from 
one  (free  space)  and  that 

|tn  - 1|  «1,  (5.2) 


where  b is  an  effective  radius;  that  the  phase  shifts  be  small.  As  a 
result  of  this  latter  requirements,  the  particle  may  be  broken  up  into 
volume  elements  (subparticles)  which  interact  only  in  terms  of  the 
phase  of  the  light  scattered  from  each  volume  element.  The  scattering 
Intensities  may  be  calculated  by  integrating  over  the  volume  of  the 
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particle  without  regard  for  any  further  interaction  of  the  internal 
electric  and  magnetic  fields. 

When  the  particles  are  large, 

a>>X,  (5.3) 

geometric  optics  is,  of  course,  valid. 

Thus  the  problem  of  light  scattering  by  particles  of  size 
approximating  the  wavelength  remains  to  be  solved.  The  modern  digital 
computer  has  permitted  several  approximate  techniques  for  the  solution 
of  this  problem  to  be  developed.  The  remainder  of  this  chapter  will  be 
concerned  with  a brief  review  of  some  of  these  techniques. 

One  of  the  earliest  of  these  techniques  was  that  developed  by 
Yeh  (Yeh,  1964).  This  technique  makes  use  of  perturbation  theory, 
expanding  the  particle  radius  (as  a function  of  angles),  the  normal 
to  the  surface,  and  the  electric  and  magnetic  fields  in  perturbation 
series.  The  boundary  conditions  on  the  fields  are  applied  and  the 
resulting  expressions  are  separated  in  orders  of  perturbation. 

Yeh's  calculations  are  limited  to  only  slightly  nonspherical 
particles  that  could  be  accurately  treated  with  first  order  pertur- 
bation theory  (preserving  the  cylindrical  symmetry  found  in  Mie 
theory).  Little  additional  work  appears  to  have  been  done  using  this 
technique,  perhaps  due  to  the  burden  of  explicitly  developing  the 
analytic  formulae  of  the  perturbation  series. 

It  is  obvious  that  for  slightly  irregular  particles  whose 
radius  differs  from  a constant  by  much  less  than  the  wavelength, 
some  type  of  perturbation  solution  is  available.  For  more  Irregular 
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particles,  it  is  necessary  that  the  Rayleigh  Hypothesis  be  Introduced. 
Because  of  its  importance,  Millar's  statement  of  the  Rayleigh  Hypothesis 
for  a two-dimensional  scalar  field  will  be  paraphrased  here  (Millar, 


1973).  The  extension  to  the  three-dimensional  vector  field  is 
straightforward . 

"Consider  a two-dimensional  configuration  (Figure  (5.1)) 
with  origin  inside  the  scatterer  S.  The  domain  exterior  to 
S is  denoted  by  D.  Let  C be  the  circle  with  center  0 
inscribed  in  S and  let  C'  be  the  corresponding  circum- 
scribed circle.  A field  is  Incident  on  S.  Then  at  all 
points  outside  C',  the  scattered  field  may  be  expressed  as 

u8  AmHf)(l)(r)exp(<m*) 

(where  Hm^(r)  ■ Hankel  function  of  the  first  kind  of 
order  m,  and  the  A^,  may  be  determined  on  C'.") 

Millar  investigated  the  validity  of  the  hypothesis  for  both  the 
bounded  particle  problem  and  the  irregular  surface  problem.  His 
conclusions  may  be  briefly  summarized  by  noting  that  boundary  value 
type  scattering  problems  may  be  validly  solved  if  the  particle  has  no 
sharp  corners  and  is  sufficiently  smooth  so  that  light  scattered  from 
one  surface  region  of  the  particle  is  not  rescattered  by  another 
surface  region  of  the  particle  (no  multiple  scattering).  All  of  the 
techniques  described  in  this  chapter  and  the  next  assume  that  the 
Rayleigh  Hypothesis  is  valid.  Any  assessment  of  the  validity  of  the 
Rayleigh  Hypothesis  is  difficult.  While  the  limits  of  its  validity  have 
been  quantified  for  certain  special  Irregular,  Infinite  surfaces,  its 
validity  for  irregular  particles  is  usually  inferred  from  comparison  of 
approximate  calculations  with  exact  calculations  such  as  the  Mie  solution. 
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One  technique  for  irregular  particle  scattering,  due  to  Reilly, 
is  a direct  extension  of  Millar's  statement  of  Rayleigh's  Hypothesis 
(Reilly,  1973).  For  simplicity,  the  description  here  will  be  limited  to 
the  scalar  wave  case.  The  actual  vector  wave  case  of  interest  follows 
directly  in  a like  manner,  but  with  additional  complexity  that  obscures 
the  methodology. 

The  scalar  Helmholtz  equation  has  the  form 

+ k^m^  ■ 0.  (5.4) 

The  solution  ^for  the  incident  field  is  known  and  the  form  of  the 

s s 

solution  $ for  the  scattered  field  is  known.  The  solution  # is 

assumed  to  be  valid  on  and  outside  the  circumscribed  sphere  around  the 
particle.  The  refractive  index  m in  the  region  outside  the  circum- 
scribed sphere  is  assumed  to  be  constant  with  the  value  m^. 

The  refractive  index  m inside  the  circumscribed  sphere  is  not 
constant  and  has  the  form 

m - m^  inside  the  particle  (5.5) 


» m^  outside  the  particle. 


(Actually  m may  vary  in  any  manner  inside  the  circumscribed  sphere, 

thus  Incorporating  refractive  index  gradients  and  voids  or  impurities 

w 

in  addition  to  the  nonspherical  shape.)  The  interior  solution  $ is 
assumed  to  have  the  form 

*W  - R(r)G(r , 6,  <J>),  (5.6) 
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where:  R ■ an  unknown  function,  and 

G » a test  function  that  Is  usually  chosen  to  be  a weighted 
sum  of  orthogonal  polynomials  so  that  R may  be  determined  without 
recourse  to  Galerkin's  method. 

The  interior  Helmholtz  equation  has  the  form 


72(RG)  + k2m2RG  - 0, 


(5.7) 


which  may  be  expanded  as 


72R  + 2VG-7R  + RV2G  + k2m2RG  - 0. 

% 


(5.8) 


Since  R Is  purely  radial,  Eq.  (5.8)  may  be  simplified  ps 


+ 2^-  Is-  + RV2G  + k2m2G  - 0. 
3r2  3r  3r 


(5.9) 


Equation  (5.9)  may  be  converted  to  a radial  (one-dimensional)  differ- 
ential equation  by  multiplying  by  G*  , the  complex  conjugate  of  G and 
integrating  over  angles. 


2 

^-|/G*Gdn  + 2||/G*|^dn+R/G*V2Gdn  + k2R/m2G*Gdft  ■-  0, 
3r  r 


(5.10) 


where  dfl  indicates  integration  over  Air.  By  defining  the  integrals 


A(r)  = /G*-^dG//G*Gdn 


B(r)  = /G*V2Gd0//G*Gdfi 


C(r)  - /m2G*GdO//G*GdfJ  , 


(5.11) 


(5.12) 


(5.13) 


t 

i $ , 

i &i 


Equation  (5.10)  may  be  reduced  to 


-MJ  + 2A(r)~  + (B(r)  + k2C(r)]R  - 0 (5.14) 

ar*- 

which  may  be  solved,  usually  numerically,  for  R(r).  Once  R(r)  has 
been  calculated,  the  scattered  field  4>  may  be  calculated  by  matching 
the  fields  at  the  circumscribed  spherical  boundary  as  shown  in  Figure 
(5.1). 

While  this  technique  is  both  elegant  and  conceptually  simple, 
it  does  not  rigorously  satisfy  the  boundary  conditions  of  electro- 
dynamics. Perhaps  for  this  reason,  this  technique  has  not  enjoyed 
extensive  use  to  date. 

Another  technique  that  has  been  demonstrated  by  Millar  (Millar, 
1973)  and  others  (Davies,  1973)  is  the  least-squares  technique  that 
is  used  to  approximately  satisfy  the  boundary  conditions  on  the 
electric  and  magnetic  fields.  The  technique  is  employed  in  the 
following  manner;  the  electric  and  magnetic  fields  are  assumed  to 
have  some  functional  form  in  terms  of  series  of  N linearly  independent 
functions  and  N expansion  coefficients  for  each  of  the  fields; 
incident,  interior,  and  scattered.  Of  course,  for  a conducting 
particle,  the  interior  fields  are  zero,  and  for  either  a conducting 
or  dielectric  particle,  the  incident  fields  are  known.  The  radius  of 
the  particle  is  assumed  to  have  some  functional  form  in  terms  of  the 
angles, 

r - f (6,  *)  (5.15) 


; 


so  that  che  surface  normal  vector  may  be  calculated 
- V[r  - f (0 , $)] 


(5.16) 


and  the  boundary  conditions  written. 


rx( E 1 + Es  - Ew)  0, 

% % % ,v  r“f(0,<i>) 


(5.17) 


nx(Hi  + Hs  _ Hw)  _/a  . - 0. 

n.  n*  'V  r = f ( 0 , d>  ) 


(5.18) 


* To  this  point,  this  technique  does  not  greatly  differ  from 
either  the  Mie  solution,  or  the  technique  to  be  presented  in  the  next 
chapter.  The  difference  will  arise  as  a result  of  the  way  in  which 
the  boundary  conditions,  Eqs.  (5.17)  and  (5.18)  are  "satisfied".  To 
illustrate  the  distinction,  consider  how  the  boundary  conditions  are 
satisfied  in  the  Mie  solution  by  requiring  that  the  equivalent  of 
Eqs.  (5.17)  and  (5.18)  be  satisfied  exactly  at  all  points  on  the 
surface  of  the  particle.  This  type  of  "satisfaction"  is  known  as 
collocative  and  is  characterized  by  no  error  being  introduced  at  the 
points  where  the  boundary  conditions  are  satisfied.  In  the  Mie  case, 
this  collocation  is  performed  analytically  so  that  no  error  is 
introduced  at  all  points  on  the  surface  of  the  particle.  Collocation 
may,  of  course,  be  used  numerically  (called  the  boundary  matching 
method)  but  this  is  done  at  a finite  number  of  points  and  error  may 
be  introduced  at  those  points  where  the  boundary  conditions  are  not 
explicitly  satisfied. 
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Least-squares  boundary  condition  "satisfaction"  introduces  an 
error  in  satisfying  the  boundary  conditions.  On  the  surface  of  the 
particle,  a number  of  points  in  the  form  of  angle  pairs,  {0^,  i.  = 

1..I,  are  selected.  Often,  these  points  are  uniformly  spaced  over 
the  surface  of  the  particle,  although  the  spacing  may  be  adjusted  to  be 
less  in  regions  of  greater  curvature,  and  greater  in  regions  of  less 
curvature.  The  accuracy  of  the  solution  is  influenced  by  this  spacing 
and  it  is  a matter  of  obvious  note  that  the  solution  is  not  unique, 
varying  with  the  placement  and  number  of  points  and  the  number  of 
functions  included  in  the  field  representation. 

These  angle  pairs  are  used  to  form  radii  { r - } and  normal 
vectors  {ji-}.  The  errors  and  c ^ are  defined  by 


c •.  - [n.*(E1+Es-EW)  ]\ 

-tl  ''v  'v,  % r~r  • * 


- ( v <»iH9-tt”>lt.r.]2. 


(5.19) 


(5  20) 


and  the  sum 


- L ♦ Ci22> 


(5.21) 


is  required  to  be  minimal  with  respect  to  the  expansion  coefficients. 

This  requirements  of  minimal  sum  square  of  errors  may  be  used 
to  formulate  the  solution  in  matrix  algebra  (Draper  and  Smith,  1967). 
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/-fl(V  V ♦i>\ 

' -fi(r2,  e2,  ^2)  1 


(5.26) 


er  *xJ. 


The  solution  may  be  found  by  the  matrix  equation 


B - (A+A)-1A+C 


(5.27) 


where  A+  ■ transpose  of  A,  and 

(A+A)  * - inverse  of  A+A  . 

Solution  of  Eq.  (5.27)  does  introduce  an  error  source  into  the 
satisfaction  of  the  boundary  conditions  compared  to  the  analytical 
collocative.  However,  the  numerical  collocative  satisfaction  would 
be  to  solve  Eq.  (5.23)  for  the  expansion  coefficients  when  I ■ N. 
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I l may  be  seen  that  the  least-squares  technique  admits  more  points 
than  the  numerical  collocative  technique  and  despite  the  inherent 
error  in  the  technique,  less  total  error  may  result,  expecially  if 
the  solution  Is  constrained  by  the  size  of  the  matrix  that  can  be 
inverted . 

Additionally  A+A  is  a symmetric  matrix  whereas  A(I»N)  is 
not  and  A+A  may  therefore  be  inverted  more  easily.  This  has  been  a 
crucial  factor  for  many  computations  in  the  past. 

In  summary,  the  least-squares  technique  does  not  exactly 
satisfy  the  boundary  conditions,  but  it  may  be  more  accurate  in  terms 
of  overall  error  than  a numerical  collocation  satisfaction,  and  it 
does  offer  some  computational  advantages.  Nonetheless,  some  additional 
error  is  introduced  because  the  boundary  conditions  are  satisfied 
numerically  rather  than  analytically. 

Another  technique  Is  the  use  of  an  integral  equation  rather 
than  a differential  equation  formalism.  Although  the  use  of  this  type 
of  formalism  is  not  new  (Shifrin,  1964),  several  recent  treatments 
deserve  special  attention  since  not  only  are  the  boundary  conditions 
accounted  for,  but  numerical  results  are  presented  (Eyres  and  Nelson, 
1976;  Uzonoglu  and  Holt,  1977;  Waterman,  1971).  Although  these  three 
formalisms  are  similar  in  development,  their  applicability  differs. 

The  formalisms  of  F.yres  and  Nelson  and  of  Uzonoglu  and  Holt  are  based 
on  the  scalar  rather  than  the  vector  wave  equation.  Alternately, 
the  formalism  of  Waterman  is  completely  general,  but  was  applied  to 
conducting  particles  only. 
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In  these  three  efforts,  the  electric  (and  magnetic)  field  is 
represented  as  a series  of  orthogonal  polynomials  and  expansion 

coefficients  that  permit  the  integral  equation  to  be  reduced  to  a 

set  of  linear  algebraic  equations  that  may  be  solved  for  the  expansion 

coefficients. 

An  interesting  modification  to  Mle  theory  to  account  for  the 
nonspher leal  character  has  recently  been  advanced  (Chylek,  Grams,  and 
Pinntck,  1976).  They  suggest  that  the  most  notable  difference  between 
spherical  and  nonspherlcal  polydisperse  phase  functions  is  the  presence 
of  glories  in  the  spherical  aerosols.  Because  the  glory  may  be 
associated  with  the  sharp  resonance  in  the  c^  and  the  d^.  in  the  region 
C^na/X  they  suggest  that  the  difference  between  the  spherical  and  the 
nonspherlcal  phase  functions  is  the  presence  (absence!  of  these 
resonances.  By  replacing  the  c^.  and  the  d^.  Mie  calculated  values  with 
the  modulo  values  in  these  regions,  agreement  with  experimental  data 
superior  to  that  enjoyed  by  Mie  calculations  has  been  demonstrated. 

Chylek  et  al.  continue  to  explain  that  the  glory  is  due  to 
surface  waves  and  that  in  an  irregular  particle,  these  surface  waves 
would  disappear.  While  the  agreement  has  been  noteworthy,  the 
validity  of  this  theory  has  been  questioned  (Kerker,  1977)  and  will  be 
addressed  in  Chapter  VIII. 

Before  concluding  this  chapter,  some  discussion  of  the  natuie 
of  the  experimental  data  is  in  order  (Holland  and  Gagne,  1970;  Plnnick, 
Carroll  and  Hoffman,  197b).  These  data  are  limited,  and  in  many  cases 
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no  supplemental  data  on  refractive  indices,  size  and/or  shape  distri- 
butions, or  polarization  are  available.  The  following  features  are 
present  in  many  cases;  the  phase  function  is  broadened  about  the  forward 
peak  relative  to  the  Mie  theory,  the  backscatter  decreases,  and  the 
amount  of  structure  decreases  but  some  structure  is  still  present. 

It  should  be  emphasized  that  even  for  data  reported  for  a monodisperse 
aerosol,  averaging  over  particle  orientation  has  already  been  conducted 
during  the  measurement  process.  As  a result  theoretical  calculations 
should  be  averaged  over  particle  orientations. 


CHAPTER  VI 

NONSPHERICAL  PARTICLE  SCATTERING  FORMALISM 

While  Che  Mle  theory  formalism  developed  in  Chapter  III  is 
applicable  to  spherical  particles  only,  the  boundary  conditions  for  £ 
and  used  to  solve  the  problem  are  general  and  applicable  to  any  shape 
particle  providing  the  Rayleigh  Hypothesis  holds.  This  hypothesis, 
reviewed  in  Chapter  V,  is  valid  if  the  mathematical  representations  of 
the  scattered  fields  accurately  describe  the  scattered  field  outside  and 
on  the  boundary  of  the  particle.  Millar  has  shown  the  validity  of  the 
hypothesis  for  conducting  particles  when  the  scattered  field  is  contin- 
uous across  the  surface  and  all  singularities  of  the  scattered  field  are 
Inside  the  particle.  For  the  dielectric  particles  to  be  considered  here, 
it  is  adequate  to  require  that  the  interior  and  scattered  fields  are  con- 
vergent and  that  the  surface  is  sufficiently  simple  that  no  secondary 
(or  higher  order)  scattering  from  one  region  of  the  particle  surface  to 
another  region  occur. 

This  discussion  will  be  concerned  only  with  nonspherical  par- 
ticles that  are  cylindrically  symmetric  with  respect  to  the  z axis,  and 
incident  fields  will  propagate  only  along  the  z axis.  These  restrictions 
allow  retention  of  the  Hertz  vector  solutions  derived  in  Chapter  III  for 
the  spherical  problem  and  in  addition,  are  forced  by  the  computational 
requirements  of  the  problem.  This  computational  restriction  is  described 
in  Appendix (I).  If  the  particle  considered  does  not  exhibit  this  symmetry. 
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or  the  incident  fields  propagate  in  some  other  direction,  higher  order 
angular  functions  must  be  included  in  the  solution  with  a proportional 
increase  in  the  number  of  expansion  coefficients  to  be  determined. 

The  surface  of  the  nonspherical  particle  is  defined  by  its 

radius 

a(0)  ■ £SnPn(cos0).  (6.1) 

The  function  describing  the  surface  is  then 


h(r,0)  ■ r-a(0)  “ 0 

- r-ES  P (0) 
n n 


(6.2) 


) 

I 


It  may  be  noted  that  the  "normalization"  or  unitization  of  the  surface 
normal  is  unnecessary  for  this  calculation. 

By  introducing  the  notation 


- S1  U“  - f 

- (f 


(6.4) 


(6.5) 
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where  the  superscripts  i,  s,  and  w indicating  the  incident,  scattered, 
and  interior  ("within")  fields  were  Introduced  in  Chapter  III,  the 
boundary  conditions  are 

* °»  (6-6) 


and 


qkaJJ  - o. 


(6.7) 


In  determinant  form.  Equations  (6.6)  and  (6.7)  may,  by  use  of  Equation 
(6.3),  be  written  as 


AE 


9 - ilPACil  AE 

30  a 9 


AE. 


and 


0. 


(6.8) 


36 

0 


AH 


6 -|J^1  AH. 


AH. 


- 0. 


(6.9) 
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Equations  (6.8)  and  (6.9)  may  be  written  in  component  form  as 


<<«£>r  - - - o 


(6.10) 


(l^Aj(J)r  “ - AH.  - 0 


30 


(6.11) 


<*“80  - - ^ " 0 


(6.12) 


^Se  " ' AH*  " 0 


(6.13) 


<r£J&  ■ »e,  + Uas“)  ‘^-o 


(6.14) 


■ H + " °- 


(6.15) 


Obviously,  Equations  (6.10)  and  (6.12),  and  Equations  (6.11)  and  (6.13) 
are  redundant,  so  that  the  boundary  conditions  reduce  to 


AL  - 0 


- 0 


(6.16) 

(6.17) 


A Eg  +1^1  AE_  - 0 


30 


(6.18) 


A^  + AH„  - 0, 


30 


(6.19) 


' 
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When  the  regular  Mie  problem  for  spherical  particles  is  solved, 
the  orthogonality  of  sin(^)  »cos(<f>)  .and  the  P^'s  is  used  to  solve  for  the 
coefficients  c^  and  d^.  Equations  (3.13)  - (3.16).  This  technique 
cannot  be  used  for  the  nonspherical  particle  problem  because  the  radial 
functions  and  Equations  (2.45) and  (2.47),  are  now  also  functions 

of  angle  0.  The  orthogonality  of  the  <j>  functions  is  still  valid  for  the 
nonspherical  particles  considered  in  this  chapter,  but  the  field  com- 
ponents will  be  seen  to  be  already  resolved  in  $,  so  that  no  advantage 
is  derived  from  this  orthogonality.  By  combining  Equations  (2.49)  - 
(2.54)  and  (3.3)  - (3.8),  the  field  of  components  for  the  three  waves 
may  be  written  as 


Er  “ + Ck.ir ) ]P^(cos0)c°s (d>) 


(6.20) 


(kr)-1£  >£ 


pJ(cos6)  3pJ(cos6) 

*Vkr)n^ — + ^<kr>as 


cos(«tl)  (6.21) 


i l 


3P^(cos0)  , Po(cO80) 

-^Vkr>ae **(kr)sin(0)  ■ 


sin(t)  (6.22) 


H* 

r 


^(kr)  + i^(kr) 


P^(cos0) sin(i) 


(6.23) 


I 


( p];(cose)  3pJ(cos9) 

Hq  “ (kr)'  lY^(kr)Bln(6')—  4-  ^(kr)~ 


sin(<p)  (6.24) 


( 3P»(cos  ) 

Hi’  (kr)"  jTc[-^(kr)30 


P»(CO80) 

^(kr)sin(0)  "J 


cos(i(i)  (6.25) 


Er  “ ^(^(kmr)  + l*'£<kmr)]p£(cos0)cos(<,)) 


(6.26) 


Eq  = (kmr) 


-1 


3P^(cos0)  P^i(cos0) 


COs(<f>) 


(6.27) 


rW 


(kmr)  ll 


l 


P.(cos0) 

-a^(kmr)gln(e)  - -t  b^(kmr)  3Q 


3P^(cos0) 


sin(<|>) 


(6.28) 


H”  - ]mhl  ^(kmr)  + ^(kmr)  P£(cos0)sin((J>) 


(6.29) 


(kr)"1! 
I 


pj(cos0)  3P^(cos0) 

^(kmr)sin(0T~t  b£^(kmr)Ii 


8ln(<|)) 


(6.30) 
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, 3P»(cos0)  P^(CO80) 

(kr)'1!  c^(kmr)^ + V*(k,,r).in(eT"  c°8^)(6*31) 


Er  " p£(kr)  +PJ  (kr>)p£<cos0)cos(<,,) 

, f 1t  9P1 (cos0)  1 Pj(cos0) 

Eq  - (kr)  l C^P£  (kr)gQ  + ^p£(kr)sin(0) 


(6.32) 


cos(4>)  (6.33) 


, pJ(cos0) 
(kr)  l -c ^ (kr)sinQ 


. 3P»(cos0) 

id^  P ^(kr)-^ sin(<t>) 


(6.34) 


Hr  ’ Id£fp^kr)  + pt"(kr)  P£(cos0)sin(4') 

D ' 


(6.35) 


? n (cos0 ) , , 3pi(cos0) 


Pft(CO80)  ,1  3P/ 

- (kr)-1!  ictp\  <kr)^ +itPt  (kr)3Q j8in<*> 


(6.36) 


, 3Pp(cos0)  ,,  Pj(cos0) 

(kr)"1!  ^(kr)^ + dtpt  (kr)8ln(e)™  coS«»(6.37) 

el 
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M 

M 


: 


— ee • r^.i 


where  the  prime  on  the  radial  functions  indicates  differentiation 


with  respect  to  total  argument,  and  ^ 


While  the  boundary  conditions,  Eqs.  (6.16)  - (6.19),  may  not  be 
solved  simply  because  of  the  angular  dependence  of  the  and  the 
functions,  these  equations  may  be  rewritten  as 


E”  - e!  - E* 


<P  <P 


(6.38 


* > 


(6.39 


Ee  - E?  + - E°)  • ES  * f^E‘ 


(6.40 


Hw  _ H,  + 


E"  - E* 


ul  . 3lna(6) ci 

H0  + 90  Er 


(6.41 


and  each  side  multiplied  by  P^,  and  integrated  with  respect  to  cos(0). 
As  an  example,  Eq.  (6.38)  becomes 


jpj,(cos0)  (E*  - E® 


dcos6 


jpj, 


(cosOJE^dcos©. 


(6.42] 


If  each  series  in  Eqs.  (6.20)-  (6.37)  is  truncated  at  t - L,  then  there 
exist  L a^,  b^,  c^,  and  d^;  4L  unknowns.  There  are  4 boundary  con- 
ditions, Eqs.  (6.38)  - (6.41),  and  integration  of  each  with  P^,  yields  L 
equations  for  a total  of  4L.  Thus  there  are  4L  equations  in  4L 
unknowns.  The  a^,,  b^,  c^,  and  d^  may  be  solved  for  by  using  some 
coefficient  of  equation  algorithm  such  as  a Gauss-Jordan  reduction 
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rc?  ■nirr- 


(Carnahan,  Luther,  and  Wilkes,  1969),  as  will  be  used  in  the  calcula- 
tions described  in  the  next  two  chapters.  It  may  be  noted  that  the  a^ 
and  the  b^  are  unwanted  quantities  as  far  as  the  desired  result  is 
concerned,  but  their  calculation  is  dictated  by  the  structure  of  the 
formalism. 

By  defining  the  integrals 


+1 


rr 


U'  ,1 


^(ka)  + ^(ka)  -^P^.dcose 


aina-1-1 


+1 


(6.43) 


kW,2 

^(ka) 

pip1, 

Ftl 

* yl 

ka 

sin(0)' 

■\p^(ka) 

x 

kU',3 

myl, 

ka 

---  D 

30  l 

’^ka) 

3PJ  „i 

kU\  4 

* yl 

ka 

30  Tl 

dcosQ 


,dcos6 


kU',5  ~ yl 


lf^'(ka)  pJpJ, 


ka  sin(0) 


dcos9 


6 ” + 'l'jj(kma) 


pJp^H^co80 


f*£(kma)  pJfJ- 

C’.7  " Ji 


UL, 7 Jkma  sin(0) 


\lCO80 


(6.44) 

(6.45) 


(6.46) 


(6.47) 


(6.49) 
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ka\  8 


i^(kma)  SP^ 
30  V 


kma 


dcos 


A££\9 


^(kma)  3P^  ^ 

30  IVdcos 


YT.10  “ J 


^(kma)  pJp^, 


kma  sin(0) 


-dcos0 


' U\ll 


| pj  (ka)  + p*"(ka)]pjpj,  H^cosO 


'll’,  12 


1 11 
c P^(ka)  P^P^, 


ka  sin(0) 


dcos0 


'll’,  13 


,-p1'  (ka)  3P»  - 

"3 r pr  dcos0 


ka 


xe;i4 


fp1»(ka)  3P»  - 

Ik5 30%*dcos6 


'll 


1 ’ 11 

rpj  (ka)  pJpJ, 

’,15  = jka  sin(0)  dcose> 


(6.50) 


(6.51) 


(6.52) 


(6.53) 


(6.54) 


(6.55) 


(6.56) 


(6.57) 


where  the  0 dependence  has  not  been  explicitly  shown  for  the  P^'s, 
the  P* 's,  and  a.  The  4L  equations  formed  form  Eqs.  (6.20)  - (6.37) 
and  (6.38)  - (6.41)  may  be  written  as 


PI 

f : 


u 
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^aZAZZ\  10  + *-bZAlt',  9 ' ClAtt\  15  ' 14>  = 


I 4 + Ate,,s) 


£ 


(6.58) 


^(-t-a^mA^,^-  b^mA^1()  + 'tc£A££'jl4  + 


| (‘^££',4  ‘ A^',5) 


(6.59) 


l(~atAlt\8  ' aZAU\6  ' 'tb£A££,,7  + CZAZZ\  1 3 + ClAtt\  1 1 


+ ^A££'12)  = lt~*AJU\2  " A^',3  " A^'^iC6-60) 


K- 


'ca£inAte',7  “b£mA££',8'b£mAt£',6  + ^CZAZJZ\  12  + d£At£',13 


= I("Atc',3  + Ate’,i  “ (6,61) 


Equations  (6.58)  - (6.61)  may  be  written  in  matrix  form  as 


B G = F 


(6.62) 


where  B is  a 4L  by  4L  matrix  given  by 
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- .... 


(6.63) 


where  each  entry  is  an  L by  L array  with  £ variation  along  the  columns 
and  variation  along  the  rows.  The  matrix  6 is  a 4L  by  1 (one)  matrix 
given  by 


(6.64) 


where  each  entry  is  an  L by  1 matrix,  and  F is  a 4L  by  1 matrix  given  by 

/ 4 + kU',S)  \ 

* kU\  \ 

F = I (6.65) 

J!(-Atfo,  . - -tA »»,  - - A ) I 

V y ■ . ' 7 

' tku\  2 ' AtP,  3/ 


where  each  entry  is  an  L (in  f')  by  1 matrix.  The  coding  of  l:.qs. 
(6.63)  - (6.65)  and  their  solution  is  described  in  Appendix  (1),  along 
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fk 

U\\0 

‘‘Atp,  9 

~kLP,  15 

’^tP,  14  \ 

'tmAU’,9 

"^££,10 

**££',  14 

Ate',is  \ 

'(Ate,.6+Ate,,8) 

7 

Atp.n*A£P,i3 

^££',12  1 ^ 

rt"lA(r,  7 

'm(A£P,6M££,,8) 

7r,i2 

Atp,ll+App,l  3/ 

with  a listing  of  the  code.  A demonstration  that  Eqs.  (6.61)  - 
(6.65)  reduce  to  Mie  theory  for  spherical  particles  is  given  in 
Appendix  (II).  Some  special  considerations  involved  in  calculating 
the  radial  functions  ^ andp|  are  discussed  in  Appendix  (III). 

After  the  c^,  and  the  d^,  have  been  calculated,  it  is  necessary 
to  form  the  and  the  3^,  used  in  the  cross  section  and  phase  function 
calculations,  Eqs.  (3.41),  (3.42),  (3.57),  and  (3.61).  These  quantities 
were  defined  in  Eqs.  (3.37)  and  (3.39)  as 


Ol  ({1*1) 


(6.66) 


and 


1-0 


t* 1 lit* 1) 


'W 


(6.67) 


It  may  be  noted  that  the  formulae  for  calculating  the  cross  sections 
and  phase  functions  are  the  same  as  those  developed  in  Chapter  111 
for  spherical  particles. 


CHAPTER  VII 


NONSPHERICAL  PARTICLE  CALCULATIONS : I 

A computer  code  was  generated  to  perform  the  calculations 
described  In  the  previous  chapter.  This  code  Is  described  In  Appendix 
(I).  Additionally,  some  special  aspects  of  the  numerical  aspects  are 
also  described,  in  Appendix  (II). 

Any  computer  code  or  algorithm  must  be  validated.  This  vali- 
dation was  approached  in  this  Investigation  by  a three  part  effort. 

The  first  part,  described  in  Appendix  (III),  demonstrates  analytically 
the  reduction  of  the  formalism  developed  in  the  previous  chapter  to  the 
Mie  solution  when  the  particle  is  spherical. 

The  second  part  of  the  validation  was  to  demonstrate  the 
numerical  equivalence  of  the  code  to  a Mie  code  for  a spherical 
particle.  This  demonstration  was  achieved  by  exercising  the  non- 
spherical  code  and  a standard  Mie  code  (Blattner,  1972)  for  the  same 
single  particles.  Approximately  fifty  cases  were  considered.  The 
c^,  d^, cross  sections,  and  phase  functions  were  compared.  In  all  cases, 
the  agreement  between  the  calculated  quantities  was  within  the  round- 
off error  of  the  machine.  This  error  was  therefore  significant  only 
beyond  the  fourteenth  figure. 

While  these  two  validation  checks  are  heartening,  they  offer 
no  real  assurance  that  the  formalism  or  the  code  is  valid  for  non- 
spherical  particles,  only  that  it  has  the  proper  limiting  behavior.  To 
continue  the  validation,  the  code  was  exercised  to  provide  calculations 
for  comparison  with  exact  calculations  for  prolate  spheroids  (Asano 
and  Yamamoto,  1975). 
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The  size  parameter  in  this  exact  theory  is 


b2  P 

1 - (7.1) 

a 


where:  a * semimajor  axis,  and 

b - semiminor  axis  of  the  ellipsoid. 

In  these  calculations,  a/b  ■ 2 in  all  cases,  and  values  of  c ■ 1,  3,  5,  7 


are  used.  The  radius  a(0)  and  its  derivative 


3a(9) 


are  calculated  both 


analytically  and  by  expansion  of  a (9)  in  Legendre  polynomials  as  a 
further  check  on  code  accuracy.  The  calculated  phase  functions  are 
shown  in  Figure  (7.1). 

Comparison  of  the  approximate  calculations  using  this  formalism 

and  the  exact  calculations  evidenced  minor  differences  in  the  fine 

structure  of  the  phase  functions.  The  positions  of  the  relative  maxima 

and  minima  and  the  relative  magnitudes  of  the  curves  agreed  very  well, 

however.  Although  it  is  difficult  to  make  an  exact  estimate  of  the 

accuracy  of  the  calculation,  an  approximate  comparison  is  made  by 

comparing  calculated  values  with  exact  values  at  10°  increments.  This 

comparison  gives  an  overall  error  bound  of  ± 2X  with  an  average  error 
_2 

of  8 x 10  %.  Because  some  of  this  error  may  be  attributed  to  the 

limited  number  of  expansion  coefficients  calculated,  to  computer  noise, 
and  to  the  truncation  error  of  the  Gauss-Jordan  integration,  this  error 
bound  seems  to  be  reasonable. 

The  overall  conclusion  of  this  validation  exercise  is  that  the 
agreement  with  limiting  cases  evidenced  by  the  code  indicates  that  it 
is  useful  despite  the  Rayleigh  Hypothesis. 
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SCATTERING  ANGLE 


FiRiire  (7. 


1)  Calculated  Phase  Functions  of  Prolate  Spheroids  Used  to 
Validate  Code  by  Comparison  with  Exact  Calculations 
(Asnao  and  Yamamoto,  1975).  (Solid  curve  ia  parallel 
polarization,  dotted  la  perpendicular  polarization.) 
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Having  validated  the  code,  attention  was  next  turned  to  the 
cross  sections  and  phase  functions  of  single,  nonspherlcal  particles. 
Individual  particles  of  shape  function 
x(6)  - 2na(0)/A 

- xo[l*O.lPn(0)],  (7.2) 

where  x - 21Ta/A,  the  Mie  parameter,  for  x ■ 1,  3,  5,  7,  9 and 
o o 

n * 2,  3,  4,  5 are  considered.  A refractive  Index  of  m ■ 1.5  + Oi 
Is  used  In  this  and  all  subsequent  calculations. 

The  calculations  are  compared  with  Mie  calculations  for  the 
same  xq.  Comparison  of  the  cross  sections  demonstrate  values  both 
greater  than  and  less  than  the  equivalent  Mie  cross  sections.  A 
similar  variation  in  the  forward  scatter  and  backscatter  is  also 
observed . 

The  phase  functions  are  also  compared  and  found  to  be  similar 

but  not  identical  to  the  Mie  phase  functions.  While  it  is  difficult 

to  draw  concise  trends  of  difference  between  the  phase  functions,  it 

is  noted  that  the  differences  increase  with  increasing  n and  decrease 

with  increasing  xfl.  These  trends  may  be  seen  in  Figures  (7.2)  - (7.9). 

In  figures  (7.2)  - (7.6),  the  phase  functions  for  x(0)  ■ xq 

(1-0. 1P^  (0)),  xo“l,  3,  5,  7,  9 are  shown  to  demonstrate  the  variation  of 

difference  of  phase  functions  with  xQ.  The  decreasing  difference  for 

xQ-l,  3,  5 is  evident,  although  the  trend  appears  to  reverse  for  xQ  “ 7 

and  9.  This  reversal  is  an  indication  of  the  inexactitude  of  the  trend. 

Figures  (7.4)  and  (7.7)  - (7.9)  show  the  increase  of  difference 

in  the  phase  functions  with  increasing  n.  In  this  case,  the  shape 

function  x(0)“  x (1-0. IP  (0)),  n ■ 2,  3,  4,  5,  and  x„  ■ 5 is  used, 
on  o 
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COSINE  OF  SCATTERING  ANGLE 


Figure  (7.2)  Unpolarized  Phase  Functions  for  Shape  Functions  x(0)  • 

x and  x(0)  ■ x (1-0. IP. (0))  (Dotted), x ■ 1. 
o 0 4 o 


SCATTERING  INTENSITY 


Figure  (7.4)  Unpolarized  Phase  Functions  for  Shape  Functions  x(6) 

5. 


and  x(0)  • xo(l-0.1PA  (0))  (Dotted), xp 
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SCATTERING  INTENSITY 


Figure 
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SCATTERING  INTENSITY 


Figure  (7.7)  Unpolarized  Phase  Functions  for  Shape  Functions  x(9)  * 

x and  x(6)  ■ x (1-0. IP,  (0))  (Dotted), x ■ 5. 
o o c o 
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Additionally,  the  effect  of  varying  the  strength  of  the  defor- 
mation is  examined.  Shape  functions  of  the  form 

x(6)  - xo[l±anPn(0)]  (7.3) 

for  n = 2 and  3,  a *■  0.05,  0.01,  and  0.15  are  considered.  Two  effects 
n 

of  this  variation  of  a are  noted;  the  values  of  the  relative  maxima  and 

n 

minima  vary  by  as  much  as  a factor  of  five  Although  the  positions  of 
the  maxima  and  minima  vary  only  slightly,  and  the  amount  of  near  forward 
scatter  (scattering  angle  less  than  30°)  varies  by  as  much  as  a factor 
of  five  about  the  equivalent  Mie  curve.  This  latter  effect  would  seem 
significant  in  terms  of  the  experimental  data  reported  by  Chylek  et  al. 
since  they  demonstrate  this  frequently  observed  effect. 
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CHAPTER  VIII 


NONSPHERICAL  PARTICLE  CALCULATIONS:  II 


In  Chapter  V,  a modification  to  the  Mie  solution  to  model  non- 


spherical  particle  scattering  was  reviewed  (Chylek,  Grams  and  Pinnick 


1976).  Because  this  modification  is  so  attractive  as  a possible  tool 


in  radiation  transfer  analyses,  part  of  this  investigation  is  devoted 


to  an  estimate  of  the  validity  of  the  modification 


This  modification  is  based  on  the  correlation  of  the  resonances 


of  the  Co  and  d»  in  the  region  x=£  with  glories  in  the  phase  functions 


glories  not  being  commonly  observed  experimentally.  As  a result 


Chylek  et  al.  proposed  the  truncation  of  the  resonances  to  model  the 


effects  of  particle  irregularity.  This  modification  was  implemented 


into  a code,  and  calculations  were  presented  that  enjoyed  better  agree 


ment  with  experimental  data  than  did  standard  Mie  calculations.  This 


truncation  was  not  applied  to  all  c»  and  d»,  but  rather  only  to  those 


for  £-3,  otherwise  the  calculations  do  not  exhibit  good  agreement  with 


the  data.  This  restriction  is  based  on  an  argument  that  a certain 


minimum  particle  size  is  necessary  for  the  validity  of  the  modification 


The  starting  point  for  this  effort  is  the  calculation  of  the 


c_  and  d in  the  region  x=£.  Shape  functions  of  the  form 


The  real  and  imaginary  parts  of  the  c and  d are  compared  with  the 


c,  and  d.  calculated  using  the  Mie  solution.  In  seven  of  these  eight 


cases,  the  shape  of  the  curves  is  preserved,  although  the  resonance 


sharpens  in  the  sense  that  the  real  parts  of  the  curve  narrow  by  as 

much  as  a factor  of  two.  This  behavior  is  demonstrated  in  Figure  (8.1) 

where  the  real  part  of  c^  is  plotted  for  the  Mie  solution  and  for  two 

shape  functions.  The  imaginary  parts  of  the  curves  also  narrow.  In 

one  case,  that  of  x(0)  = x (1+0.1P_( 0)) , the  resonance  curve  is  replaced 

o J 

by  a slowly  varying  curve  whose  real  parts  have  a mean  of  approximately 


0.3,  and  whose  imaginary  parts  have  a mean  of  approximately  zero. 

Additionally,  although  the  centers  of  the  curves  shift,  apparent 
ly  without  trend,  the  real  parts  of  the  curves  always  fall  under  the 
Mie  resonance  curves.  The  extent  of  the  narrowing  or  the  curves  seems 
to  be  approximately  proportional  to  n.  On  the  basis  of  these  calcu- 
lations, therefore,  the  modification  of  the  c^  and  d^  is  not  sub- 
stantiated, but  the  numerical  results  of  Chylek  et  al.  may  be  justified 
in  view  of  the  narrowing  of  the  resonance  curves. 

As  a further  step,  calculations  are  performed  for  one  of  the 
polydisperse  aerosol  distributions  reported  by  Chylek  et  al.  (i.e.  the 
first  KC1).  Nine  phase  functions  are  calculated  using 
x(0)  - x 


x ( 0 ) = X (l±0.05Po(6)) 
o i 


x ( 0 ) = x (ltO.lPo(0)) 
o l 


x(0)  = x (1*O.15P~(0)) 
o i 


x(0)  - x (1— O.lPo(0)*O.lP. (0)) 


Mie  Parameter  x 


Figure  (8.1)  Comparison  of  real  parts  of  c_  for  x(0)  * x (1-0.1P„(6)) 

x(0)  = x (1-0. IP 0 (0) ) I ),  and  Mie  (solid). 


and  combined  linerly  with  adjustable  coefficients  calculated  using 
regression  techniques  to  fit  the  data.  Two  different  fitting  schemes 
are  used;  one  using  all  the  data,  and  another  using  only  those  data 
for  scattering  angles  greater  than  sixty  degrees.  A representative 
fit  for  each  scheme  is  compared  with  the  data  in  Figure  (8.2).  These 
fits  appear  to  be  as  good  as  those  reported  by  Chylek  et  al. 

It  must  be  noted,  however,  that  the  Chylek,  et  al  modification 
incorporates  the  effect  of  averaging  over  particle  orientation  while 
the  formalism  of  Chapter  VI  is  limited  to  only  one  particle  orien- 
tation. A complete  calculation  would  incorporate  this  averaging 
process.  Unfortunately,  this  average  is  beyond  the  scope  of  this 
investigation  because  of  the  limitation  of  available  computer  memory. 
If  the  averaged  values  of  c ^ and  were  calculated,  better  agree- 
ment with  the  behavior  suggested  by  Chylek,  et  al.  might  be  observed. 
Additionally,  improved  agreement  with  the  experimental  data  shown  in 
Figure  (8.2)  for  scattering  angles  between  100  and  170  degrees  might 
result.  It  is  in  this  region,  however,  that  the  resonance  behavior 
of  the  c^,  and  d^  is  most  Important.  While  the  agreement  between 
these  calculations  and  the  data  is  poorest  in  this  region,  the  agree- 
ment is  still  better  than  that  displayed  by  Chylek,  et  al's  calcu- 
lations because  truncating  the  resonances  in  the  c^  and  d^.  effectively 
removes  all  structure  such  as  is  found  in  the  experimental  data  in 
this  region  of  scattering  angle.  Thus,  some  resonant  behavior  may 
still  be  expected  in  the  eg  and  dp  after  averaging  over  particle 
orientation,  and  the  behavior  suggested  by  Chylek,  et  al,  should  not 
be  completely  reproduced. 


Figure  (8.2  ) 


Comparison  of  NonspherleaJL  Data  (Circles)  of  Chytek, 

Grams,  and  Pinnick  with  Figurative  Mie  (Solid)  and  Two 
Combinations  of  Calculated  Phase  Functions.  Both 
combinations  use  phase  functions  for  the  shape  functions 
x(0)  ■ xo(1-0.05P2(6)),  x(9)  ■ x (1-0. 1P2(0)) , and  x(0)  ■ 

xQ(l-O.15P2(0)) . The  curve  ( ^ was  fitted  using  all 

data,  and  the  curve  (-.-)  was  fitted  using  only  data 
for  scattering  angles  greater  than  sixty  degrees. 


l * , 
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CHAPTER  IX 


REVIEW  OF  RADIATIVE  TRANSFER  THEORY 


The  theory  of  Radiative  Transfer  (RT)  has  been  developed  to 
describe  the  transport  of  light  in  an  extinguishing  medium.  The 
general  four-dimensional  form  of  the  RT  equation  is  (Chandrasekhar, 
1960;  Pomraning,  1973) 


the  j component  of  the  Stokes  vector  (Eqs 


where 


■ position  vector, 

J^,  Jj'  ■ unit  propagation  vectors, 

a(£,  t)  - the  extinction  coefficient  at  % and  t (Eq  (4.2)) 
u)(r,  t)  ■ single  scattering  albedo  at  £ and  t (Eq.  (4.6)), 
p..^r,  t,  k*k’)  ■ scattering  phase  function  for  the  /’ th 


Stokes  vector  component  into  the  J 
Stokes  vector  component  at  J^and  t through 
an  angle, 

the  cosine  of  the  scattering  angle. 


The  phase  function  *yy  tCfc*  f»  differs  from  the  P - (^,  t. 


defined. in  Eq.  (4.4)  by  a pair  of  rotations  to  correct  for  the  change 


of  orientation  from  to  Jg'  (Chandrasekhar,  1960). 


The  RT  equation  is  most  commonly  used  in  its  unpolarized, 


time  independent  form 


a(r)oj(r)  r 

vpcf  t>  - t>  ♦ «i Jpce.  v>d  V 


(9.2) 


where:  I(^,  j^)  ■ the  unpolarized  intensity,  and 

p(£,  J^’J^')  “ phase  function  defined  by  Eq.  (4.5). 
Additionally,  a source  function  J(,£,  J^)  is  commonly  defined  as 


J<t-  <t>  -*?}»<*.  «'»<*•  V>dV. 


(9.3) 


which  allows  Eq.  (9.2)  to  be  rewritten  as 


V " Jfe>] 


(9.4) 


Until  recently,  most  of  the  work  on  RT  has  been  performed  in 


describing  stellar  or  planetary  atmospheres  (Chandrasekhar,  1960), 


This  work  was  concerned  with  plane  or  spherical  parallel  atmospheres 


with  symmetry  of  medium  and  boundary  conditions  that  allowed 


Eq.  (9.4)  to  be  reduced  (in  the  plane  parallel  case)  to 


y)  - -a(Z)[l(Z,  u)  - u>(Z)  J (Z,  y)] 


(9.5) 
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where  Z = coordinate  along  which  the  medium  varies 


U “ Jj.*Z,  and 


J(Z,  M) 


P(Z,  Jfc*Jfc')I(Z,  p’)dp'. 


Equation  (9.5)  is  commonly  referred  to  as  the  RT  equation  in  the 
literature  because  of  the  interest  it  has  enjoyed. 

While  Eq.  (9.5)  is  useful  for  describing  the  brightness  of  a 
star  or  the  amount  of  sunlight  reaching  the  surface  of  a planet,  it  is 
limited  to  media  varying  along  only  one  axis  and  uniform  incident 
illuminations.  Thus  Eq.  (9.5)  cannot  be  used  to  treat  problems  where 
the  media  varies  in  more  than  one  direction  such  as  smoke  (Friedlander , 
1977;  Greene  and  Lane,  1964)  or  clouds  (Mason,  1971).  Additionally, 
nonuniform  incident  illuminations  such  as  are  encountered  in  imaging 
and  laser  propagation  problems  cannot  be  treated. 

To  date,  calculations  performed  on  varying  media  and/or 
nonuniform  incident  illuminations  have  been  limited  to  either  few 
orders  of  scattering  (Deepak  and  Green,  1970)  or  small  angle  scatter- 
ing (Weinman  and  Shipley,  1972).  Additionally,  some  Monte  Carlo  RT 
calculations  have  been  performed  for  laser  transmission  through 
clouds  (Bird,  1974),  but  these  calculations  are  costly  in  terms  of 
computer  time  required  because  the  trajectory  of  each  photon  must 
be  followed  from  source  to  absorption  by  either  the  medium  or  the 
detector,  and  a large  number  of  photons  must  be  counted  to  achieve 


"•  ■ .LL._  .■  ipJK§.! 


u -.» ..  » 


statistical  significance.  Thus,  there  have  been  no  general  numerical 
RT  calculations  performed  to  date  from  a general  solution  of  the  RT 
equation. 


The  RT  equation,  Eq.  (9.1)  may  be  rewritten  by  noting  that  the 


,1  3 


differential  (—  ^ + J^*^)  Is  essentially  a directional  derivative 


Je 


where 


l ’ Zo  + # 


(9.6) 


t = t + a c 
o 


r^  = boundary  of  the  medium,  and 


t * time  the  light  crosses  the  boundary. 


Equation  (9.1)  may  be  rewritten  as 


t'  k)  m ~ a(S’  + 


(9.7) 


a<*>  t)ai(jr,  t)  , * 2 

P<t.  t,  V)d\ 


4tt 


in  an  unpolarized  form.  Because  most  experimental  data  provide 
only  the  information  of  a concentration  of  total  aerosol  as  a 
function  of  position  and  time  and  a constant  particle  size  distrib- 
ution, it  is  most  useful  to  adopt  the  case  where  N(^,  a,  t)  is  separ- 
able (Eqs.  (4.1)  - (4.4)),  so  that  Eq.  (9.7)  may  be  rewritten  as 
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3/*^’  *■»  a^»  c»  fO 


(9.8) 


a(£,  t)w 


4tt 


p^*^')1^*  c>  fc,)dV 


th 


If  the  scattering  events  are  Independent,  n order  scattering 
is  due  to  n consecutive  single  scatterings,  and  the  intensity  may  be 
written  as  a sum  of  contributions  for  each  order  of  scattering  (Deepak 
and  Green,  1970).  That  is. 


N 


x0e>  c»  fc) " c*  # 


(9.9) 


n*0 


th 


where  I (j^.t,  k)  is  the  intensity  at  £ at  time  t along  £ due  to  n 
order  scattering.  While  it  is  assumed  that  this  series  is  convergent 
it  may  be  noted  that  the  series  will  converge  in  the  limit  n -*■  “ if 
wl  1 (£,  t,  J^)/Ir(^,  t,  J^)  < 1.  Since  w<l,  this  assumption  seems 
reasonable  and  is  borne  out  numerically. 

Equation  (9.8)  can  be  decomposed  into  N equations  by 
substituting  Eq.  (9.0); 


3lV*»  C’  & 


(9.10) 


and 
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aZV*’  c’  V = -a(^’  t)In(*’  t*  J6) 


a(^,  t) 


)In_l(^»  c»  fc’>dV»  n>°  (9‘U) 


It  is  convenient  to  extend  the  definition  of  the  source 
function  by  combining  Eqs.  (9.3)  and  (9.9)  to  give 

Jn(*’  C’  ‘ i?  jp^V^V*’  C’ 


(9.12) 


Equation  (9.11)  may  be  rewritten  as 

IfV*’  c’  V = "a(*’  t)In(^‘  fc*  V + a(*’  c) Jn-1(^’  C’  & 


(9.13) 


The  boundaries  of  the  medium  may  be  expressed  as  the  set  of 
vectors  {r^}  and  the  initial  times  as  {tQ}.  The  incident  illuminations 
are  defined  as  I0 ( » 1 0 )•  The  solutions  to  Eqs.  (9.10)  and  (9.13) 
may  then  be  written  as 

l 

V*>  c’  V “ V**’  V JS)  exp(  ~j  to+£'/c)d£’j(9.14) 


V*>  c*  k)m  In(^o*  *=0’  V exp  ("  to+£'/c)d£’]  (9.15) 


aW’  to+r/c)Jn-lW’  to  + V/c»'k) 


exp  - 


vr/c)d£" 


v 


dV  , n>0. 


If  the  boundaries  are  nonreflecting,  t * Jj)  = 0;n>0,  and  Eq. 


(9.15)  may  be  reduced  to 

l 


Vt'  V ■ 


“W-  '0+e,cV  (9-16) 


exp 


f 

- a(r  +k£",  t +l"/c)dl" 
%o  o 


dl'. 


V 


Equations  (9.9),  (9.12),  (9.13),  and  (9.16)  constitute  the 
four  basic  equations  for  the  general  RT  calculations  to  be  performed. 
Chapter  X describes  the  development  of  a numerical  algorithm  for  the 
solution  of  the  two-dimensional  RT  equation,  and  Chapter  XI  presents 
some  calculations  performed  using  the  code  implementing  this  algorithm. 
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CHAPTER  X 


TWO-DIMENSIONAL  RADIATIVE  TRANSFER  EQUATION  ALGORITHM 


The  RT  equation  is  solved  by  introducing  a quadrature 


approximation  for  the  source  function  similar  to  that  of  Chandrasekhar 


(Chandrasekhar,  1960).  The  approximation  for  the  source  function 


where  the  W.  are  approximate  weighting  factors  chosen  such  that 


where  ev  is  the  unit  vector  in  the  z direction.  The  phase  function 


is  now  independent  of  position  and  time  because  the  medium  has  been 


assumed  to  be  separable.  Additionally,  because  the  RT  equation  of 


interest  in  this  chapter  will  be  essentially  two-dimensional,  all 


reference  to  the  functional  dependence  on  time  t will  be  dropped 


only  steady  state  conditions  will  be  admitted 


The  aerosol  medium  will  be  represented  by  a rectangular  array 


L by  M points  in  extent.  The  positions  of  the  points  in  the  array  are 
designated  by  the  indices  l and  m.  Since  the  array  is  rectangular, 


the  points  in  the  array  have  eight  first,  second,  and  third  nearest 


neighbors,  defining  eight  directions  (and  J^) . This  geometry  is 
shown  in  Figure  (10.1). 


Mmiiuni, r .mi1  wu^  iin-mii.j i.  h-m.j  ■ ip  'hp^pwi* 


» 


I 


2 

The  differential  d in  Eq.  (9.12)  may  be  written  as 
dcosO,d4>'  and  Kq.  (10. 1)  may  be  formed  as  the  product  of  two  trapezoid 
rule  Integrations  (without  derlvitive  terms) (Abramuw It z and  Stegun, 
19<>4).  One  property  of  the  trapezoid  rule  is  that  the  weighting 
factors  for  all  points  except  the  end  points  are  equal.  The  weighting 
factors  for  the  end  points  are  one-half  of  the  other  weighting  factors. 
The  integration  is  approximated  at  the  points  4''  “0,  n , and  2n. 
Since  the  integrand  of  Kq.  (9.12)  has  the  same  value  at  4*'  “ 0 and  at 
4>'  - 2ir  , the  integration  effectively  tooks  only  at  the  points  41'  ■ 0 
and  4>'  - n with  equal  weighting  factors  at  both  points.  The  cost)' 
Integration  is  approximated  at  the  points  cosO*  - ±1,  i 0.5  and  0, 
but  the  end  points,  cost)'  ■ tl  are  counted  twice  because  of  the  41' 
approximation.  As  a result,  the  weighting  factors  are  the  same  for 

all  < and  are  defined  by  reducing  Eq.  (10.2)  to 


I 


l 

w)  P(®3*&t')  “ (10.3) 

t - 1 

where  W is  the  integration  weighting  factor. 

Equation  (10.1)  allows  the  solution  of  the  time-independent 
RT  equation  to  he  reduced  to  an  eight-beam  formalism,  so' that  only 
those  intensities  needed  for  the  source  functions  are  calculated. 
Additionally,  the  approximation  that  the  source  functions  and 
intensities  are  constant  in  the  rectangular  region  about  each  point 
in  the  arrav  is  made. 
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where  Qy^  - Wp  and  ^(>)  is  the  position  in  the  array.  It  is 

conceptually  convenient  to  introduce  the  notation 


In^'^,  " In,fm,t’ 


(10.5) 


(10.6) 


and 


a<W  ■ °V 


(10.7) 


The  incident  illuminations  are  specified  at  all  points  on  the 
boundary  of  the  array  (t  ■ 0 and  ( - L for  all  ffl,  and  m » 0 and 
m - M for  all  l)  for  all  k interior  to  that  £,m  point.  The  I » 

4,  O | , J 

ure  calculated  rrom  these  boundary  conditions  by  a point-to-point, 
along  each  path  Integration  of  Eq.  (9.14).  ’’’he  Jq  ^ • are  cal- 


culated from  the  I » ■ using  Eq.  (10.4).  Following  this,  the 

o,  t m,j 


I,  ~ are  calculated  from  the  J • using  Eq.(9.16)  and  the 

l.cm ,J  o,  m,j 

J,  „ are  calculated  from  the  I,  « . . The  intensities  I « ■ 

1 , cm . j 1,  vn,j  n.twt.j 


and  the  source  functions  J -for  n>l  are  calculated  in  a similar 

n.OM,  j 


manner.  The  total  Intensity  at  any  point  in  any  direction  may  be 
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calculated  from  Eq.  (9.9).  For  example,  the  forward  intensity  along 
the  boundary  ( ■ L is 


«.v 


l I 

n-0 


0) 


n,  Lm,  1 


(10.8) 


and  the  backseat tered  intensity  along  the  boundary  ( ■ 1 is 


N 

I(*l  ,*e3}  “ ^Q^.lm.S 


M> 


These  two  quantities  are  calculated  to  the  order  of  scattering  N for 
which  the  Center  Line  of  Sight  (CLOS  - the  line  between  the  source- 
detector  system,  being  along  varying  £ for  m - M/2)  forward  intensity 
varied  with  order  of  scattering  by  less  than  IX.  That  is, 


N 

J. 


N-l 


n»o 


fM 

2 


wn  - l I fu  ,wn<.01 1 I „ <an. 

• 1 S-o  "•  1 n-o  1 


(10.10) 


Another  quantity  calculated  is  the  modulation  contrast  defined  by 


" I^m,^>m.-xI(^m.^)»in 


(10.11) 


where  I(*£m,  ^t)max  and  I(*£m,  ^)min  are  the  maximum  and  mlnlmum 
Intensities  along  the  path  perpendicular  to  Specifically,  the 

contrast  for  the  forward  intensity  along  m for  l m L is  calculated 
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since  this  would  be  the  contrast  seen  by  a detector  looking  along  the 
CLOS. 

Available  computer  memory  would  only  permit  this  two- 

dimensional  algorithm  to  be  coded  and  then  by  storing  only  the  J » , 

n § w«  i j 

the  J , . , • and  the  forward  scattered  and  backscattered  intensities, 

n-1,  cm,  J 

Eqs.  (10.8)  and  (10.9).  The  run  time  of  the  code  was  minimized  by 

introducing  an  eight  beam  interlace  integration  scheme  for  the 

intensities  and  source  functions. 

The  eight  beam  interlace  Integration  scheme  arises  from  noting 

that  the  quadrature  approximation  to  the  source  function,  Eq.  (10.4) 

may  be  calculated  sequentially  by  incrementing  for  each  1^  ^ • 

as  it  is  calculated  rather  than  waiting  until  all  eight  1^  ^ j 

have  been  calculated  and  then  summing  them.  The  chief  advantage  of 

this  scheme  is  that  it  allows  the  calculation  of  the  source  functions 

to  be  shifted  from  a consideration  of  each  point  in  the  array  on  a 

point  by  point  basis  to  a consideration  of  each  point  in  the  array  on 

a path  by  path  basis,  collecting  points  as  they  occur  on  each  path 

to  increment  the  source  function  at  that  point.  The  operation  of 

this  scheme  is  shown  in  Figure  (10.2). 

The  advantage  of  Eqs.  (9.9),  (9.12),  (9.14)  and  (9.16)  in 

coding  is  that  they  allow  the  intensity  for  each  order  of  scattering 

to  be  calculated  from  the  source  function  for  the  previous  order 

of  scattering.  The  code  may  be  structured  so  that  the  Iq  ^ ■ are 

calculated  from  the  boundary  conditions,  the  Jq  ^ • are  calculated 

from  the  I » and  in  general,  the  I . • are  calculated  from 

o.  On, j n,  un,j 


§ j 

: •*.  A 
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the  J , „ ; and  the  J ; from  the  I . ...  The  calculation  of 

n-l.&n,)  n,  lm,j  n,  V”,  J 

the  1 ; has  been  described  by  Eq.  (10.4).  The  equivalent 

tl  | vn  t J 

equations  for  the  calculation  of  the  I « • from  the  incident 

illuminations  and  the  ^ j from  the  Jn_^  fa  j be  derived 

from  Eqs.  (9.13)  and  (9.16),  respectively.  To  describe  this,  consider 
the  specific  example  of  the  first  two  points  along  the  direction 
jL  ■ 1.  These  two  points  are  t ■ 1,  and  l ■ 2,  ■ 1.  Equation  (9.13) 

may  be  rewritten  as 


Io,  li,  i " Io,  oi,  iexp  ^“a  nAZ/2-* 


(10.12) 


and  since  there  are  no  reflecting  boundaries,  and  the  source  function 
is  constant  in  the  rectangular  region  about  each  point  in  the  array, 
Eq.  (9.16)  may  be  rewritten  as 

AZ/2 

11.  1 ' “llVl,  11.  1 f Kl<“'2  * 4Z',ldM' 


(10.13) 


which  may  be  integrated  to  yield 


rn,  11,  1 “ Vl.  11,  1 Cl  - CXP  <-allAZ/2>3 


(10.14) 


The  upper  limit  of  integration  AZ/2  occurs  since  this  is  the  distance 
from  the  boundary  of  the  array  to  the  point  l ■ 1,  w ■ 1 along  the 
direction  i ■ 1.  The  intensities  Iq  ^ ^ and  In  j each 


I 


contribute  one  term  Q^Iq  1 aru*  Qyj^n,  11,  1 to  C^e  ca^cu^at^on 

of  J_  ..  ; and  J , respectively.  The  Intensities  at  the 

U|  1 1 i J n i iifj 

point  ( ■ 2,  im  - 1 may  also  be  calculated  form  Eqs.  (9.13)  and 


(9.1b)  us 


AZ  3AZ/2 

V 21.  1 ■ V 01.  1 Ud4Z'  - | 421d“' 

0 AZ 


-In  ,1  lexP  “ailAZ  " °21AZ/2 


(10.15) 


lo,  li,  r 


AZ  ar 

21,  1 ‘ “llVl.  11.  1 H'  I'  I “l 


(10.16) 


3AZ/2 


3AZ/2  3AZ/2 


-L  iSZ"  ♦a21Jn_i.  21  ! Jddz'exp-  a21dS2"j, 

(_  A*w  A7« 


which  may  be  integrated  to  yield 


V 2i.  i - Vi.  n.  1(i-«p(-»ll4»)-P^2iaz/2) 


+Vi.  a.  iHV1™ 


(10.17) 
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since  the  point  l ■ 2,  m » 1 Is  3AZ/2  from  the  boundary,  and  the 
source  function  changes  value  Is  going  from  the  first  rectangular 
region  about  / - 1,  m - 1 to  the  second  rectangular  region  about 
f - 2,  m ■ 1.  Equation  (10.15)  may  be  rewritten  as 


Xn,  21,  1 " rn.  11,  1 eXp  *"(ail  + a2i)AZ'/2l 


(10.18) 


by  recognizing  that  Iq  ^ ^ is  equal  to  the  intensity  Iq  ^ ^ 
degraded  by  extinction  between  the  points  l ■ 1,  si  ■ 1,  and  L m 2, 
m ■ 1,  while  Eq.  (10.17)  may  be  rewritten  as 

In,  21,  1 " Jn-1,  11,  1 t1-exP(-a11Az^2HexP[-(a11+oi21)AZ/2] 


(10.19) 


+Jn-1  11  [l-expt-cij^  ^AZ/2)  ] exp  (— ot0  AZ/2 ) 


+Jn-1,  21,  1l1-exPt“a21AZ/2^ 


The  first  term  of  the  rtght  hand  side  of  Eq.  (10.19)  may  be  recognized 

as  the  intensity  l . . . degraded  by  extinction  between  the  points 

n t Al|  1 

l - 1,  in  - 1 and  f » 2,  m - 1.  Equation  (10.19)  may  then  be  rewritten 


1 - I . exp[-(a  -W  )AZ/2] 

n,  21,  l n,  11,  1 11  21 
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+Jn-1  11,  ill-exP(-a11AZ/2) )exp(-a21AZ/2) 

{ 


+Jn-1  21  1 U-exp(^2iAZ/2)  1 . (10.20) 


Equations  (10.18)  and  (10.20)  are  the  specific  forms  for  this  example 

of  the  general  numerical  equations  for  the  intensity  at  a point  in  the 

array  along  a path  in  terms  of  the  intensity  at  the  previous  point  in 

the  array  along  the  path  and  the  previous  order  of  scattering  source 

functions  at  the  two  points  on  the  path  along  the  path,  noting  that 

I » “0  for  n^0  for  all  (,  m,  /.  It  mav  also  be  noted  that 

n,  tin, j 

Eq.  (10.20)  may  be  derived  from  Eq.  (9.15)  if  ^ and  Jfc)  “ 

I ..  , . It  is  the*,  combination  of  Eqs.  (10.4),  (10.12),  (10.14), 

n,  11,  1 

(10.15)  and  (10.20)  that  allows  the  order  of  scattering  path-bv-path 
calculation  to  be  performed. 

The  derivation  of  these  equations  incorporated  several 
approximations.  Two  approximations  were  used  in  deriving  Eq.  (10.4); 
that  the  extinguishing  medium  may  be  represented  as  planar,  and  that, 
within  that  plane,  the  radiation  transport  may  be  represented  by  eight 
beams.  It  is  expected  that  computational  errors  will  increase  when  the 
Incident  illumination  varies  greatly  outside  the  plane,  or  when  eight 
beams  are  not  enough  to  represent  accurately  either  the  incident  i 
illumination  or  the  scattering  phase  function.  Equations  (10.12), 
(10.14),  (10.15),  and  (10.20)  incorporate  one  major  approximation, 
that  the  source  function  does  not  vary  appreciably  about  each  point 
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in  the  rectangular  array  representing  the  planar  medium.  The  iterative 
form  of  the  solution  and  the  close  computational  connection  between 
Eqs . (10.4),  (10.12),  (10.14),  (10.15),  and  (10.20)  makes  an  exact 
estimate  of  error  difficult. 

As  an  estimator  of  worse  case  error  to  be  expected  from  the 
combination  of  these  approximations,  calculations  were  performed 
using  this  two-dimensional  algorithm  (a  listing  of  the  code  is  given 
in  Appendix  (IV))  and  a three-dimensional  Monte  Carlo  trajectory 
tracing  code  (House  and  Avery,  1969)  for  a Rayleigh  scattering 
medium  (Chandrasekhar,  1960)  with  an  albedo  of  0.5.  The  medium  was 
assumed  to  be  uniform  (ot^,  * constant).  The  Rayleigh  phase  function 

was  chosen  for  this  comparison  because  it  represents  the  worse  case 
of  any  situation  of  interest  with  regard  to  nonforward  scattering, 
since  polydisperse  aerosol  phase  functions  do  not  generally  exhibit 
glories  or  rainbows  (Van  de  Hulst,  1957)  Deirmendj ian,  1964).  The 
incident  illumination  was  a gaussian  profile  beam  held  constant  with 
respect  to  the  third  dimension  for  the  Monte  Carlo  calculations. 

The  forward  CLOS  intensities  at  five  optical  depths  are  given  in  Table 
(10.1).  The  precent  differences  were  calculated  relative  to  the 
Monte  Carlo  results.  It  may  be  seen  that  while  there  are  considerable 
differences  between  the  two  calculations  (maximum  3.3%),  the  variation 
in  the  difference  and  the  degree  of  approximations  extant  in  the 
two-dimensional  algorithm  indicate  that  the  agreement  in  this  case 
may  be  considered  to  be  quire  good.  Better  agreement  might  be 
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expected  for  most  polydisperse  aerosol  phase  functions,  but  worse 
agreement  for  incident  illuminations  with  more  variation  outside  th« 
plane  when  the  planar  medium  approximation  becomes  more  inaccurate. 


CHAPTER  XI 


RESULTS  OF 

TWO-DIMENSIONAL  RADIATIVE  TRANSFER  CALCULATIONS 

The  calculations  presented  in  this  chapter  are  for  a uniform 

Deirmendjian  C.3  fog.  Because  the  fog  is  uniform,  the  extinction 

coefficient  is  constant  within  the  fog  which  encompasses  the 

entire  rectangular  array.  In  terms  of  Eq.  (4.9),  D = 5.5556, 

3 

B = 0.3333,  6 = 3,  Y =8,  and  the  concentration  is  100/cm  (Deirmendjian 
1964).  Three  wavelengths  were  considered;  1.06  ym,  3.0ym,  and  10.6  ym, 
representing  the  near  ir,  the  mid  ir,  and  the  far  ir.  Two  of  these 
wavelengths  correspond  to  actual  laser  wavelengths  (1.06  ym  and  10.6  ym) 
The  other  wavelength  was  chosen  close  to  the  3-5  ym  window  of  the 
atmosphere,  but  just  outside  the  window  where  the  water  absorption  would 
be  considerable  (Selby  and  McClatchey,  1976).  Additionally,  the 
refractive  index  of  water  is  well  known  at  these  wavelengths  (Centeno, 
1944;  Pointer  and  Dechambenoy,  1966;  Curcio  and  Petty,  1954). 

The  extinction  coefficients,  albedoes,  and  phase  functions 
were  calculated  using  an  existing  polydisperse  Mie  scattering  code 
adapted  to  calculate  the  0 • • (Blattner,  1972).  A listing  of  this  code 
is  given  in  Appendix  (IV).  The  wavelengths,  refractive  indices, 
calculated  extinction  coefficients,  albedoes,  and  Qy^'s  for  / " 1-5 
are  given  in  Table  (11.1).  The  other  Qy^'s  are  not  included  in 
Table  (11.1),  since  the  symmetry  of  the  unpolarized  phase  function 

requires  that  Q21  = Q81'  q31  * Q71’  Q41  " Q61* 

Two  sources  were  considered  in  these  calculations  — a uniform 

intensity  beam  and  a Gaussian  profile  beam.  The  uniform  intensity 


beam  was  used  to  investigate  the  transmission  of  contrast.  This  beam 
was  1.9  km  wide  against  a black  background  5.2  km  wide  including  the 
beam.  The  geometry  of  these  calculations,  showing  the  uniform  intensity 
beam,  is  given  in  Figure  (11.1).  As  defined  by  Eq.  (10.11),  this  beam 
lias  an  Inherent  contrast  of  one.  Both  beams  were  treated  as  plane 
waves  traveling  along  the  i ■ 1 direction  for  the  purpose  of  establish- 
ing boundary  conditions.  The  other  boundary  conditions  were  taken  to 
be  zero.  For  both  beams,  the  incident  illumination  at  the  CLOS  was 
one  so  that  intensities  along  the  CLOS  could  be  directly  interpreted  as 
transmissions. 

The  depth  of  the  fog  was  varied  in  0.5  km  Increments  from  0.5 
to  5 km.  The  size  of  the  array  was  M - 31  and  L - 50,  limited  only  by 
the  available  computer  memory.  The  scale  of  Az  is  determined  by  the 
depth  of  fog  to  be  considered,  the  scale  of  Ax  is  determined  by  the 
scale  of  Az,  and  the  size  of  the  beam  is  determined  by  the  number 
of  Ax  needed  to  accurately  represent  it.  Thus  the  beam  diameter  is 
directly  constrained  by  the  depth  of  fog  considered. 

The  forward  and  backscattered  intensities  are  calculated  using 
Equations  (10.8)  and  (10.9)  for  both  beams  at  the  three  wavelengths 
and  ten  fog  depths  described.  The  forward  intensities  along  the  CLOS 
are  shown  in  Figures  (11.2)  and  (11.3)  for  the  uniform  profile  and 
gaussian  profile  beams,  respectively.  Examination  of  these  intensities 
reveals  that  the  intensities  fall  off  with  fog  depth  z as  approximately 
exp  (-aqz)  where  q is  a diffusion  exponent  (always  t 1)  similar  to  that 
found  in  traditional  one-dimensional  plane-parallel  RT  (Chandrasekhar, 
1960;  Kattawar  and  Plass,  1967).  The  intensity  at  3.0pm  falls  off 
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Figure  (11.3  ) Forward  CLOS  Intensity  of  Gaussian  Beam. 
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fastest  because  the  extinction  coefficient  is  largest  at  that  wave- 
length due  to  the  water  absorption.  While  these  q's  are  not  the  same 
numerically  as  those  encountered  in  one-dimensional  RT  because  of  the 
nonuniformity  of  the  boundary  conditions,  this  behavior  of  the  inten- 
sity along  the  CLOS  is  not  unexpected  because  of  the  symmetry  of  the 
boundary  conditions.  The  intensity  of  the  uniform  beam  falls  off 
faster  than  the  intensity  of  the  gaussian  beam  for  each  wavelength,  and 
the  magnitude  of  the  difference  between  the  intensities  for  each 
wavelength  decreases  with  increasing  wavelength.  This  difference  in 
intensities  is  due  to  both  scattering  and  initial  beam  profile.  The 
total  orders  of  scattering  N required  to  satisfy  the  convergence 
condition.  Equation  (10.10)  is  dependent  on  the  depth  of  fog  and  the 
albedo,  which  is  wavelength  dependent.  As  the  albedo  or  the  fog  depth 
increases,  N increases.  Experience  with  the  two-dimensional 
code  indicates  that  N is  more  strongly  influenced  by  albedo  value 
than  by  fog  depth.  Further,  the  value  of  q decreases  with  increasing 
u)  and/or  N.  This  behavior  of  q has  been  noted  in  plane-parallel 
one-dimensional  RT  (Kattawar  and  Plass,  1976). 

The  difference  in  transmission  between  the  two  beams  is  also 
due  to  the  form  of  the  initial  beam  profile.  The  uniform  beam  inten- 
sity falls  off  less  rapidly  adjacent  to  the  CLOS  than  does  the  gaussian 
beam  incident  intensity.  This  condition  affects  the  scattering  of 
light  away  from  the  CLOS.  The  total  scattering  out  from  the  CLOS  is 
less  for  the  uniform  beam  than  for  the  gaussian  beam  for  low  orders 
of  scattering  since  the  intensity  scattered  out  of  the  CLOS  is  largely 
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tile  beam,  intensity  that  scatters  out  into  the  background  is  not 
replaced,  largely  because  there  is  little  intensity  to  be  scattered 
from  the  background  into  the  beam.  Some  intensity  is  scattered  into 
the  edge  of  the  beam  from  the  interior  of  the  beam,  but  there  is  a 
net  scattering  out  from  the  beam  into  the  background.  As  higher  orders 
of  scattering  are  considered,  more  intensity  scatters  into  the  back- 
ground than  out  of  it. 

This  net  scattering  out  of  the  CLOS  also  occurs  for  the 
gaussian  beam,  but  at  a lesser  rate  since  there  are  no  regions  of  zero 
initial  intensity  for  the  gaussian  beam.  If  the  gaussian  beam  were 
truncated  to  the  size  of  the  uniform  beam,  the  intensity  of  the 
truncated  gaussian  beam  would  fall  off  even  faster  than  the  intensity 
of  the  uniform  beam. 

The  forward  intensity  is  strongly  influenced  by  the  total 
orders  of  scattering  N required  to  satisfy  Equation  (10.10).  As  pre- 
viously stated,  the  value  of  N required  is  largely  determined  by  the 
albedo.  In  practice,  it  is  found  that  the  relative  value  of  N 
required  goes  at  least  as  the  square  of  the  albedo.  Since  the  value 
of  the  albedos  used  in  this  study  decrease  with  increasing  wavelength. 

N decreases  with  increasing  wavelength.  Thus  shorter  wavelengths 
evidence  higher  order  intensities  that  fall  off  faster  for  the  uniform 
beam  than  for  the  gaussian  beam. 
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Comparison  of  the  backscattered  intensities.  Figures  (11.4) 
and  (11.5),  respectively,  at  each  wavelength  shows  that  the  uniform 
beam  intensities  are  greater  than  the  gaussian  beam  intensities. 

Further,  the  magnitudes  of  the  intensities  decrease  with  increasing 
wavelength  as  do  the  differences  between  the  intensities.  This  is 
because  the  backscattered  intensity  is  largely  due  to  low  orders  of 
scattering.  This  fact  is  demonstrated  by  the  shape  of  the  curves, 
being  proportional  to  1 - exp(-2aq ' z) , where  q'  is  again  a diffusion 
exponent.  This  curve  may  be  calculated  analytically  for  first  order 

A 

scattering,  using  Equation  (9.15)  for  n - 1 and  ^ - -e^.  While  q' 
could  not  be  calculated  as  exactly  as  q for  each  wavelength-beam  case 
because  of  the  difficulty  of  accurately  calculating  the  proportionality 
constant  for  each  curve  (essentially  a total  source  term,)  it  was  found 
that  q'  < q in  all  cases  by  an  amount  consistent  with  the  errors 
involved  in  the  IX  criterion  of  Equation  (10.10).  This  behavior  is 
consistent  with  the  contention  that  backscattered  intensity  is  largely 
due  to  low  orders  of  scattering  since  the  1%  condition  applies  to  the 
forward  intensity  and  for  equal  numbers  of  forward  and  backscattered 
intensity  terms  contributing  to  the  total  intensities,  the  backscattered 
Intensity  should  be  more  accurate  and  q'  < q. 

The  increase  in  backscatter  with  decreasing  wavelength  may  be 
seen  by  considering  that  the  albedo  decreases  with  increasing  wave- 
length, and  the  product  of  albedo,  extinction  coefficient,  and 
decreases  with  increasing  wavelength.  This  latter  quantity  is  the 
backscatter  coefficient  that  determines  the  amount  of  first  order 
backscatter  back  along  the  CLOS.  Thus  the  backscatter  coefficient 
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Figure  (11.5  ) Backscatter  CLOS  Intensity  of  Gaussian  Beam. 
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determines  the  first  order  baekscatter,  and,  in  this  case,  the  relative 
wavelength  order  of  the  curves,  while  the  albedo  determines  N and 
therefore  i}'. 

The  relative  difference  between  the  backscattered  intensities 
at  each  wavelength  is  due  to  the  same  effects  that  caused  the  differ- 
ences in  the  forward  intensities;  scattering  and  initial  beam  profile. 
The  tirst  order  backscattered  intensity  along  the  CLOS  is  the  same 
for  both  beams.  The  second  and  higher  order  backscattered  intensities 
are  due  to  intensity  from  outside  the  CLOS  entering  the  CLOS  to  be 
scattered  back  along  the  CLOS.  Since  these  entering  intensities  fall 
off  as  exp(-al),  the  immediate  region  to  the  CLOS  is  most  important  in 
contributing  to  the  backscattered  intensity.  For  the  uniform  beam, 
this  region  changes  less  in  intensity  for  low  order  scattering,  and 
more  for  high  order  scattering  than  it  does  for  the  gaussian  beam. 

Thus,  since  the  backseat tered  intensity  is  largely  determined  by  low 
orders  of  scattering,  the  uniform  beam  should  have  more  backscattered 
intensity  than  does  the  gaussian  beam. 

The  scattering  of  intensity  out  of  the  uniform  beam  is 
demonstrated  in  the  contrast  of  the  uniform  beam  forward  intensity. 

The  contrast  for  the  three  wavelengths,  calculated  using  Equation 
(10.11),  is  shown  in  Figures  (11. o)  - (11.8),  respectively.  Since 
the  contrast  is  symmetric  about  the  CLOS,  only  one-half  of  the  beam 
contrast  is  shown.  Figure  (11. b),  the  contrast  of  the  uniform  beam 
at  1.06;m  clearly  shows  the  decrease  of  the  contrast  of  the  beam,  the 
increase  of  the  contrast  of  the  background,  and  the  deformation  of  the 
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beam.  These  three  features  are  most  pronounced  in  this  figure,  and 
may  be  seen  in  Figures  (11.7)  and  (11.8)  to  decrease  in  magnitude 
with  increasing  wavelength.  Figure  (11.7),  the  contrast  of  the 
uniform  beam  at  3.0pm,  shows  this  decrease  in  magnitude,  although 
the  decrease  >of  the  beam  contrast  and  the  increase  in  the  back- 
ground contrast  are  still  evident  in  the  beam  contrast  curve  furthest 
from  the  CLOS  in  crossrange  and  in  the  background  contrast  curve 
nearest  the  CLOS,  respectively.  The  deformation  of  the  step  func- 
tion shape  of  the  Initial  profile  is  much  less  evident  at  3.0pm  than 
it  was  at  1.06pm.  At  10.6pm,  Figure  (11.8),  increase  of  the  back- 
ground contrast  is  barely  dtscernable,  although  some  decrease  in 
beam  contrast  and  deformation  of  beam  profile  may  be  seen  at  extreme 
fog  depth.  , 

The  gaussian  beam  profiles  also  demonstrate  the  scattering 
of  intensity  outward,  although  to  the  lesser  degree  expected.  For 

this  reason,  one-half  of  the  beam  profile  at  1.06pm  for  the  gaussian 
beam  is  presented  in  Figure  (11.9).  This  profile  was  «a<lculafced  by 

dtviding  the  Intensity  at  each  fog  depth  and  crossrange  by  the  CLOS 
intensity  at  that  fog  depth.  The  deformation  of  the  beam  profile  to 
a wider  gaussian  is  evident  in  the  increase  of  the  relative  inten- 
sity curves  away  from  the  CLOS  with  increasing  fog  depth.  The 
deformation  may  be  seen  to  be  much  less  for  this  gaussian  beam  than 
in  the  comparable  uniform  beam.  Figure  (11.6).  The  beam  profiles 
at  3.0pm  and  10.6pm  are  not  shown  since  they  differ  only  qualitatively 
from  Figure  (11.9).  Figure  (11.10)  shows  the  increase  of 
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the  beam  diameter  with  fog  depth  at  the  three  wavelengths.  The 
relative  effects  of  albedo  value  may  be  clearly  seen. 
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CHAPTER  XII 


SUMMARY 


This  dissertation  has  described  the  results  of  an  investigation 


of  radiation  transfer  through  aerosols.  The  first  part  of  this  inves 


tigation  was  concerned  with  single  particle  scattering.  The  scattering 


of  light  by  a single  spherical  particle,  the  Mie  solution,  was  reviewed 


by  developing  the  Hertz  vector  solution  for  the  Maxwell  Equations  in 


Chapter  II  and  reproducing  the  Mie  solution  using  the  Hertz  vectors 


in  Chapter  III.  The  extension  of  the  Mie  solution  to  describe  the 


single  scattering  of  light  by  a collection  of  independent  particles 


the  usual  case  with  aerosols,  was  reviewed  as  developed  by  Deirmendjian 


in  a format  compatible  with  performing  radiative  transfer  (RT) 


calculations  (Deirmendjian,  1964).  This  description  of  aerosol  single 


scattering  was  given  in  Chapter  IV 


The  Mie  solution  is  valid  only  for  spherical  dielectric 


particles  so  that  it  may  treat  most  liquid  aerosols  and  those  solid 


aerosols  that  are  spherical.  Large  liquid  aerosol  particles  that  have 


been  deformed  by  gravity  and  aerodynamic  drag,  and  most  solid  (irregu- 


lar) aerosol  particles  cannot  be  treated.  While  dielectric  particles 


of  regular  geometric  shape  may  be  treated  by  solving  the  Mie  problem  in 


other  coordinate  systems,  no  generally  accepted  solution  for  irregular 


particles  is  available.  Much  effort  has  been  directed  towards  this 


problem,  however,  and  several  solutions  have  been  developed  to  treat 


various  aspects  of  the  irregular  particle  problem.  Several  of  these 


The  central  theme  of  the  Irregular  particle  scattering  problem 
is  the  Rayleigh  Hypothesis  which  limits  the  degree  of  irregularity 
of  the  particle  if  the  scattered  fields  encorporate  only  outgoing 
spherical  waves.  Thus  solutions  that  assume  the  Rayleigh  Hypothesis 
are  not  valid  for  particles  that  are  sufficiently  irregular  for  light 
to  be  scattered  from  one  part  of  the  particle  to  another  and  then  out. 
All  of  the  solutions  described  in  Chapter  V assume  the  Rayleigh 
Hypothesis. 

The  previous  efforts  described  in  Chapter  V satisfy  the  bound- 
ary conditions  of  electrodynamics  approximately.  The  implication  of 
these  approximations  is  unclear.  Recently,  an  integral  equation 
solution  has  been  advanced  that  accounts  for  the  boundary  conditions 
correctly  and  solves  the  vector  Helmholtz  equation  rather  than  the 
scalar  (Waterman,  1971).  Unfortunately,  only  conducting  particles 
have  been  treated  numerically  with  this  method  although  it  is  valid 
for  dielectric  particles. 

Additionally,  Chapter  V reviews  the  Mie  solution  modification 
proposed  by  Chylek  et  al.  (Chylek,  Grams  and  Pinnick,  1976).  This 
modification  is  based  on  the  observation  that  glories  are  not 
commonly  observed  in  the  experimentally  measured  phase  functions  of 
polydisperse  irregular  particles.  The  mathematical  source  of  glories 
in  the  Mie  solution  is  the  resonances  in  the  scattered  field 
expansion  coefficients.  Chylek  et  al.  propose  that  these  resonances 
be  truncated  to  model  the  effects  of  particle  irregularity  on 
scattering.  While  the  validity  of  this  modification  has  been 
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questioned,  Chylek  et  al.  have  demonstrated  better  agreement  with 
experimental  data  for  their  calculated  phase  functions  that  for  Mie 
phase  functions. 
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Recognizing  the  limitations  of  previous  efforts  to  treat  the 
problem  of  irregular  particle  scattering,  the  Mie  theory  was  extended 
to  cylindrically  symmetric  irregular  dielectric  particles.  This 
extension  is  presented  in  Chapter  VI.  The  limitation  of  this  extension 
to  cylindrically  symmetric  particles  was  not  prompted  by  an  deficiency 
in  the  basic  theory,  but  rather  by  the  constraints  of  available 
computer  memory.  The  only  theoretical  limitation  on  this  solution  is, 
in  common  with  all  other  solutions  described  in  this  dissertation,  the 
assumption  of  the  Rayleigh  Hypothesis,  which  permits  only  a small 
degree  of  particle  irregularity. 

A computer  code  was  generated  to  calculate  the  scattering 
properties  of  these  particles.  This  code  is  listed  in  Appendix  (I). 
Analytic  and  numerical  calculations  were  performed  to  demonstrate  that 
the  theory  and  code  would  reproduce  the  Mie  solution  when  spherical 
particles  were  treated.  These  calculations  are  described  in  Appendix 
(III)  and  Chapter  VII,  respectively.  The  Mie  solution  has  also  been 
extended  to  nonspherical  particles  of  regular  geometric  shape.  Notable 
among  these  efforts  is  the  solution  for  prolate  and  oblate  ellipsoids 
(Asano  and  Yamamoto,  1975).  Approximate  calculations  using  the 
solution  of  Chapter  VI  were  performed  for  prolate  ellipsoids  arid 
compared  with  the  exact  calculations.  Despite  the  Rayleigh  Hypothesis, 
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the  agreement  evidenced  by  this  comparison  was  quite  good.  In  view 
of  this  agreement,  the  applicability  of  the  code  to  large  irregular 
particles  is  demonstrated.  These  calculations  are  described  in 
Chapter  VII. 

This  code  was  then  used  to  Investigate  the  scattering  proper- 
ties of  single  irregular  particles,  and  to  address  the  validity  of  the 
Mle  solution  modification  of  Chylek  et  al.  Calculations  were  performed 
to  determine  the  effect  of  particle  irregularity  on  the  resonant 
behavior  of  the  scattered  field  expansion  coefficients.  While  it  was 
found  that  the  resonances  narrowed,  the  numerical  method  of  Chylek 
et  al.  was  not  substantiated.  This  finding  did  not  explain  the  agree- 
ment of  the  Chylek  et  al.  calculations  with  the  data,  so  the  code  was 
modified  to  treat  polydisperse  aerosols  using  the  techniques  described 
in  Chapter  IV.  Calculations  were  then  performed  to  generate  poly- 
disperse phase  functions  for  nine  shapes  of  particle  using  the 
particle  size  distribution  and  refractive  index  reported  by 
Chylek  et  al.  for  one  set  of  their  data.  These  nine  phase  functions 
were  then  combined  linearly  in  various  combinations  using  regression 
techniques  to  fit  the  data.  It  was  found  that  several  combinations 
of  calculated  phase  functions  would  fit  the  data  as  well  as  the 
calculations  of  Chylek  et  al.,  and  that  the  choice  of  phase  functions 
used  in  fitting  the  data  was  not  especially  crucial.  The  insensi- 
tivity of  the  data  to  particle  shape  indicates  that  calculations 
utilizing  their  modification  may  be  adequate  for  most  RT  calculations. 
Chapter  VIII  is  devoted  to  a description  of  these  calculations. 
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The  RT  theory  uses  certain  single  scattering  quantities,  the 
extinction  coefficient,  the  albedo  of  single  scattering,  and  the  phase 
functions,  along  with  the  appropriate  boundary  conditions  to  permit 
solution  of  the  problem  of  describing  the  multiple  scattering  effects 
of  aerosols.  This  theory,  mainly  developed  by  Chandrasekhar,  had 
boon  primarily  used  to  treat  the  problems  associated  with  stellar  and 
planetary  atmospheres  where  the  symmetry  and  time- independence  of 
the  boundary  conditions  and  the  extinguishing  medium  allow  the  general 
four-dimensional  (3  space  plus  time)  RT  equation  to  be  reduced  to  a 
one-dimensional  RT  equation  (Chandrasekhar,  1960).  This  theory  and 
the  development  of  the  one-dimensional  equation  are  reviewed  in 
Chapter  IX. 

The  RT  theory  has  most  commonly  been  applied  to  practical 
problems  such  as  the  calculation  of  the  brightness  of  a star  or  the 
amount  of  light  that  reaches  the  surface  of  a planet.  A new  set  of 
practical  problems,  that  of  describing  the  propagation  of  low  energy 
laser  beams  and  images  through  the  atmosphere  with  and  without  aerosol 
aerosols,  has  developed  in  recent  years  that  cannot  be  treated  by 
one-dimensional  RT  theory.  These  problems  are  characterized  by  non- 
uniform  boundary  conditions  and/or  asymmetric  extinguishing  media. 
Previous  efforts  to  treat  these  problems  have  included  Monte  Carlo 
trajectory  tracing  codes  that  are  exceedingly  time  consuming,  small 
angle  approximation  developments  of  the  four-dimensional  RT  equation 
that  are  valid  for  aerosols  only  in  the  visible  and  near-ir  regions 
of  the  spectrum,  and  are  not  valid  at  all  for  atmospheres,  and 
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general  angle  analyses  that  are  limited  to  few  orders  of  scattering. 

No  general  analytic  or  numerical  solution  of  the  four-dimensional 
RT  equation  has  yet  been  developed. 

In  Chapter  X,  a numerical  algorithm  for  the  solution  of  the 
two-dimensional  RT  equation  with  nonuniform  boundary  conditions  and 
asvmmetric  medium  is  developed.  This  algorithm  is  suitable  for 
treating  steady  state  laser  beam  and  image  propagation  problems. 

The  extinguishing  medium  used  in  this  algorithm  may  incorporate  both 
aerosol  and  atmospheric  effects  and  the  solution  may  be  obtained  to 
an  arbitrary  order  of  scattering.  Thus,  two  limitations  of  previous 
efforts,  the  small  angle  approximation  and  the  restriction  to  few 
orders  of  scattering,  are  removed. 

A code  was  generated  implementing  this  algorithm.  A listing 
of  the  code  is  presented  in  Appendix  (IV).  While  this  code  requires 
a large  amount  of  computer  memory  to  execute,  it  is  rapidly  executable. 
This  speed  advantage  was  demonstrated  when  the  code  was  validated  by 
comparison  with  a Monte  Carlo  trajectory  tracing  code  that  was  much 
smaller  in  size  but  required  much  more  execution  time.  The  accuracy 
evidenced  in  the  comparison  of  the  results  of  the  two  codes  and 
the  speed  of  execution  of  the  two-dimensional  code  demonstrate  the 
utility  of  the  code  for  steady  state  problems. 

The  two-dimensional  code  was  then  used  to  calculate  the  effects 
of  aerosol  multiple  scattering  on  transmission,  backscatter,  and 
contrast  transmission  through  an  aerosol  medium.  Calculations  were 
performed  for  a liquid  (spherical)  aerosol,  a Deirmendjian  C.3  fog 
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of  uniform  concentration.  Two  beams,  one  of  uniform  profile  and  the 
other  of  gaussian  profile,  were  considered  at  wavelengths  in  the  near, 
mid,  and  far-ir.  These  calculations  are  presented  in  Chapter  XI. 

In  summary,  two  parts  of  the  problem  of  radiation  transfer 
through  aerosols  have  been  addressed  in  this  dissertation  investi- 
gation. The  first  part,  that  of  single  scattering  by  a single 
article,  resulted  in  an  extension  of  Mie  theory  to  nonspherical 
particles.  This  extension  is  limited  only  by  the  ubiquitous  Rayleigh 
Hypothesis  and  the  current  state-of-the-art  in  available  computer 
memory.  The  code  implementing  this  extension  offers  the  capability 
to  treat  those  solid,  irregular  aerosol  particles  that  could  not  be 
treated  by  the  Mie  solution. 

While  the  second  part  of  the  radiation  transfer  problem,  that 
of  single  scattering  by  a collection  of  particles,  was  not  addressed 
because  most  aerosols  are  comprised  of  independent  particles  and  the 
basic  work  of  Deirmendjian  is  valid,  the  third  part  of  the  problem, 
that  of  multiple  scattering,  was  addressed.  A numerical  solution 
algorithm  for  the  two-dimensional  RT  equation  was  developed  and  a 
computer  code  implementing  it  produced.  This  code  is  useful  for 
treating  steady  state  problems  of  low  energy  laser  beam  and  image 
propagation  through  aerosol  laden  atmospheres. 

This  dissertation  has  suggested  some  interesting  problems 
for  future  effort.  One  problem  would  be  the  comparison  of  the 
multiple  scattering  effects  of  liquid  (spherical)  and  solid  (non- 
spherical) aerosols.  An  example  of  this  problem  is  the  difference  in 
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propagation  through  fog  as  compared  to  dust.  This  comparison  could 
be  addressed  by  using  the  single  scattering  quantities  calculated 
with  the  nonspherical  scattering  code  in  solving  the  two-dimensional 
RT  equation. 

This  second  problem  is  an  extension  of  the  two-dimensional  RT 
solution  algorithm  to  four  dimensions.  The  two-dimensional  algorithm 
developed  during  the  dissertation  investigation  is  limited  because  it 
cannot  treat  time  dependent  problems  such  as  calculating  the  stretching 
effect  of  multiple  scattering  on  laser  pulse  propagation,  nor  problems 
in  which  the  third  spatial  dimension  is  asymmetric.  While  the  two- 
dimensional  algorithm  is  adequate  for  some  modern  theories  of  detection 
and  recognition,  the  four-dimensional  RT  equation  solution  is  needed 
to  consider  problems  with  two-dimensional  images  and  vertically 
varying  media  such  as  smoke  and  dust  clouds.  The  extension  of  the 
two-dimensional  algorithm  may  be  effected  by  returning  to  the  four- 
dimensional RT  equations  developed  in  Chapter  IX.  A numerical 
solution  algorithm  may  then  be  developed  for  these  equations  by 
using  the  quadrature  approximation  and  evaluating  the  path  integrals 
numerically  rather  than  analytically  as  is  done  in  the  two-dimensional 
algorithm. 
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APPENDIX  I 


NONSPHERICAL  PARTICLE  SCATTERING  CODE 


The  computer  code  described  in  this  appendix  was  generated  to 
implement  the  formalism  developed  in  Chapter  VI.  Calculations  performed 
using  this  code  are  presented  in  Chapters  VII  and  VIII. 

The  integrals  Eqs . (6.43)  - (6.57)  are  performed  with  a 

Cuuss-Legendre  Quadrature,  (Carnahan,  Luther,  and  Wilkes,  1969). 

Equation  (6.62)  is  solved  for  the  a^,  b^,  c^,  and  d^,  with  a Gauss- 
Jordan  elimination  routine  modified  to  treat  complex  entries  (Sellers 
and  Gibbs,  1977).  The  radial  functions  ij^(x)  and  p|,(x),  Eqs.  (2.45)  and 
(2.47),  and  their  derivatives  are  calculated  with  a special  routine 
Implementing  the  method  of  Miller  (Miller,  1950).  This  method  is 
described  in  Appendix  (II).  The  cross  sections  and  phase  functions 
are  calculated  in  the  same  manner  as  the  Mie  solution  (Deirmend j ian, 

1964;  Van  de  Hulst,  1957) . 

The  formalism  developed  in  Chapter  VI  is,  of  course,  also  valid 
for  particles  not  cylindrically  symmetric  with  respect  to  the  incident 
plane  wave.  At  this  time,  however,  an  extension  of  the  formalism  and 
computer  code  to  consider  more  general  geometries  is  not  feasible  due 
to  computer  limitations.  If  a general  particle  and/or  incident  direction 
were  to  be  considered,  Eqs.  (3.3)  - (3.8)  would  Include  terms  cos(m<t>), 
sin(m^),  and  P^(cosO),  m - 0 to  L.  The  size  of  the  array  B,  Eq.  (6.63), 
would  increase  by  terms  »n  and  rw* . Arrays  6 and  F,  Eqs.  (6.64)  and 


(6.65),  would  increase  in  size  from  4L  to  4L  . Thus  the  total  storage 

2 4 2 

required  would  increase  from  16L  + 8L  to  16L  + 8L  . For  large  value 
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of  L,  this  is  an  approximate  increase  by  a factor  of  L • The  increase 
in  computer  memory  required  for  operational  code  (as  compared  to  array 

storage)  should  be  approximately  proportional.  The  present  code  requires 

5 2 4 

2x8  words  to  initiate  compilation,  so  that  L ~8  words  and  a general 

9 

code  would  require  approximately  2x8  words  to  initiate  compilation. 

Computer  memory  of  this  extent  is  available  only  on  a very  few  machines 

so  that  the  implementation  is  not  feasible  at  this  time. 

It  should  be  noted  that  while  the  code  requires  a large  amount 

of  computer  memory  for  initial  compilation  (2  x 8"*  words),  this  does 

not  restrict  the  utility  of  the  code  as  much  as  might  be  expected. 

This  code  is  normally  executed  on  the  MICOM  CDC  6600  computer  using  the 

SCOPE  3.4.2  compiler.  This  is  a two  pass  compiler  that  performs  some 

optimization  of  the  compiled  code.  In  this  case,  the  optimization  is 

significant  as  the  computer  memory  required  for  the  compiled  code  is 
4 

only  6.5  x 8 words.  Additionally,  execution  of  the  compiled  code  is 
accomplished  in  about  15  seconds  for  a nominal  particle  with  19  expansion 
coefficients  in  each  Debye  potential,  and  Mie  parameter  (xq)  of  nine. 
Compilation  of  the  code  requires  about  60  seconds.  Thus,  the  compiled 
code  may  be  used  to  perform  calculations  with  only  moderate  demands  on 
computer  memory  and  execution  time. 

This  code  consists  of  one  driver  routine  and  five  subroutines. 
These  routines  are  listed  in  Tables  (1.1)  - (1.6).  The  operation  of 
each  routine  is  described  below. 
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NONSPH  (Table  (I. 1) ) 

This  is  the  main  driver  routine.  It  accesses  the  five  sub- 
routines to  perform  the  calculation  of  the  nonspherical  particle 
scattering.  The  version  of  the  code  shown  here  is  the  polydisperse 
aerosol  code  with  provision  for  a log  normal  particle  size  distri- 
bution. The  routine  calls  START  to  input  the  necessary  parameters  and 
load  the  arrays  with  the  points  and  weights  for  the  Gauss-Legendre 
quadrature.  The  value  of  the  radius  and  its  derivative  at  these  points 
are  calculated.  The  increments  to  the  A^,,^  integrals  are  calculated 
and  the  routine  calls  ADO  to  store  these  increments  in  the  appropriate 
places  in  the  B and  F arrays.  Once  the  B and  F arrays  have  been 
calculated,  CGAUSS  is  called  to  solve  for  the  a £,  b^,  c^,  and  d^. 
Routine  AP2  is  then  called  to  calculate  the  cross  sections  and  phase 
functions  from  the  c^  and  d^. 

START  (Table  (1.2)) 

This  subroutine  performs  two  functions.  It  first  loads  the 
values  of  the  weights  and  points  for  the  Gauss-Legendre  quadrature, 
and  then  inputs  the  necessary  parameters  for  the  calculation.  These 
parameters  include  the  refractive  index,  the  wavelength  of  light,  the 
minimum  and  maximum  particle  radii,  the  average  particle  radius  and 
standard  deviation  for  the  log  normal  particle  size  distribution,  and 
the  number  of  terms  in  the  particle  radius  expansion,  the  number  of 
expansion  terms  in  the  Debye  potentials,  and  the  number  of  particle 
radii  to  be  calculated.  The  routine  also  calculates  the  increment  of 
particle  radius. 


AP2  (Table  (1.3)) 


Subroutine  AP2  performs  the  calculation  of  the  cross  sections 
and  phase  functions  for  both  the  individual  particles  and  for  the 
polydlsperse  aerosol.  Additionally,  it  calculates  the  fraction  of 
particles  in  each  radius  using  a log-normal  particle  size  distribution. 

ADO  (Table  (1.4)) 

This  routine  receives  the  increments  of  the  A,,,, ,,  integrals 

u n 

to  form  the  B and  F matrices.  This  routine  is  called  once  for  each 
•ingle  in  the  quadrature  for  each  t and  C pair. 

CGAUSS  (Table  (1.5)) 

Subroutine  CGAUSS  is  modified  version  of  a library  Gauss- 
Jordan  elimination  routine  (Sellers  and  Gibbs,  1977).  The  only 
modifications  performed  to  it  were  to  permit  consideration  of  complex 
rather  than  only  real  numbers. 

BESSL  (Table  (1.6)) 

This  subroutine  implements  the  method  of  Miller  described  in 
Appendix  (II)  to  calculate  the  radial  functions  tpg(x)  and  ^(mx). 

The  radial  function  P^(x)  and  the  derivatives  of  the  radial  functions 
iJ’j. (x),  i^(mx),  and  P^(x)  are  calculated  using  the  method  described  in 
Appendix  (11).  This  routine  has  an  upper  limit  on  x due  to  the  use 
of  the  backwards  recursion  method  of  Miller  (Miller,  1950).  Should 
larger  values  of  x be  required,  another  version  of  BESSL  impelmenting 
both  Miller’s  method  and  the  standard  forward  recursion  method  is 


available.  This  additional  capability  is  not  generally  required 
because  values  of  x that  result  in  divergence  were  not  needed  in 
this  investigation. 


F ' 

j 

f I 


ii 


' 

I 


I 


r 


THIS  PAGE  IS  BEST  QUALITY  PRACTICAL 
yjyyi  QQfi  FURNISHED  TO  LOG  — 


»! 


I 

t 

j 

I 

I 

J 

i 


i 


i 

i 


i 

i 

i 

J 

i 

! 


TABLE  (1.1).  LISTING  OF  ROUTINE  NONSPH 


s 


10 


IS 


£0 


<*s 


30 


3S 


AO 


as 


SC 


SS 


PROGRAM  NorsPH<INP'ir.OUTPUT.TAPES*iNPUI.T*PFf>»i>UTP  r> 

COMMON/ONt/X (201 .w(?0)  .AM.alAm.nm.nn.am  IS) 

COMMON/TWO/AIBO.Hl) .Cl (15> »»! 

01 PENSION  RK  120) .Rh<?0) 

COMMON/TMRtE/AM(20)  .AJ(20>  . lAJICO)  ,l)AH(,>ll)  • 'rtn<->UI 

COMPLEX  AJ.AM.PN.OAJ.OAH.OHri 

COMMDN/FOUR/NT.NP.PMIN.OR.SO.SR.AXX  < IS) .NXX 

COMPLEX  A 

COMPLEX  Cl 

COMPLEX  AM.RKM 

COMPLEX  A 1 * ALP 

OAT  A P)/3.1a1S<»27/ 

C 

C NON-SPHE^lCAL  PART  ICLE  MIF  SCATTERING 
C ASCO  PROGRAM 

C COilEO  MY  b • M.  FOwLtR 
C VERSION  OF  1 AilGUSI  IV77 

C METMOO  OS  SOLUTION  IS  HOUNnARY  VAL'F  Of  IANGFNTIAL  t.  A Ij  M fill  'is 

C USING  COEFFICIENT  of  EQUATION  sOLUlION 

C REFERENCE  IS  HORN  AND  WOLF 

C PROJECTION  USING  LFotNORE  PL0VNOM1  .'.I.S  of  FIPS1  ASS'ClAffu  I YPf 

C COt;E  LlMlfF O TO  PARTICLES  wlTH  CYLINDRICAL  XYH"M». 

C RAl.'IuS  SURFACE  EXPANUEO  in  /FroIh  I FOtooRE  PLOYNOMiALS 
C EXTERIOR  KtGION  ASSUMEO  VACUUM 

C • .»•••<■••«  ••■••••••••« 

C X = POINTS  IN  GAUSS-LEGtNORt  QUADRA TUKt 

C W * WEIGHTS  IN  GAUSS-LFGENOPF  UiIAOhATijpE 
C RK  s RADIUS  AT  XII) 

C A IS  MATRIX  Of  COFFFICIFNf  INItGRAi  s 

C Cl  IS  SllRRtNT  CONTrIhUT  ION  TO  L.N  INltoWAL  I 
C NM  « NOMritw  OF  LEGENDRE  CONMPONf  NT IN  pmOIiiS 
C ALAM  IS  WAVELENGTH  In  MICRONS 

C 

C SUBROUTINE  TO  LOAD  VARIOUS  P4RAMF  If  rs 

CALL  START 

C CALCULATION  OF  K VALUt  »WAvtNJM~.t-l 

AA*2.*PI/ALAM 

C A I IS  SQUAnt  ROOT  OF  -1 

Al«CMPLX(0..l.) 

WRITFIS.NIIAM.A1. An 

C COi't  TO  CALCULATE  T«t  VALII.S  Of  Tf  «A;.|us  AT  >■  ACM  g-L  I N f.  Ii^l  I 1 uN 
00  Nh  NsI.nI 

AN  1 1 )*RMlN. FLOAT  In- 1 I *0R 
IF  (NXX.LT.^I GOTO  4/> 

00  A)  1 I *«? .NXX 
4|  AN  < I I I * AN ( i I *AXX  < I I) 

A?  CONTINUE 
NMsNXX 
NP  = N 

00  1 I » 1 .«?U 

C P YARlAHLtS  ARF  /f -UIH  LEG.  Nt)Rf  s 

PO  = I. 

P|-=X(II 

C T vAwIAHLEs  ARF  Off-lVAUVIS  OF  /»  IM  LFOtNIIRLi  - am  > M'ilsT  m R 

T0  = 0. 

T 1 * X 4 1 1 

C PP  VARIABLES  APf  VAN  Ot  HULST  > PI  FNCs 


I 


I 
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os 


/•i 


8* 


VS 


100 


t OS 


llo 


PP0«0. 

PPt«l. 

PICll)«AN<l>*P0*AN<2)BPl 

PH(l>«AN<2l*T| 

00  2 J»3.N^ 

AJJ»FLUAT<J) 

P2bM2.*AJJ-1.1*X1II*PI*AJJ*P0(/(A  IJ-l.l 

PO*P| 

PI*P? 

PPb ( ( 2 . B A JJ- 1 . ) • A C I I »PP 1 -A J J«PPn I / ( A JJ- 1 . ) 

TT**(l)»(PP-PPn»-«?.»AJJ-| .»•(!. -XI  I >•*< !> (APPl.TO 

PP0*PP1 

PP|*PP 

1 0 = T I 

ti*tt 

HK(I)bpK(1)*AN(J)*p2 
PmU>»*<M(1>»ANIJ>»TT 
2 CONTINUE 

UK  < I I Bp*  ( I ) >A8 
l CONTINUE 

VI  F0nM«T<SX»2E20.H> 

ZERO  OUT  InIEGPAL  AkHAy 
NNNba«nN 
NNNPbNNNM 
00  3 LIbI.NNN 
00  3 N l b l « nNNP 
< a«l1.ni)«cnpl*co..o.i 
call  EuuAUON  POUT  I Nt 

MUST  EITHLk  BE  COMPLEX  OP  00  IN  T«'<  bTAC.ES 
00  A 1*1.20 

INTEGRATION  LOOP  - pEFEkS  lu  10  PI  i.-L  OUAOPAMwF 

PKMbuK(I)«AM 

AM  IS  HFFnACTIVE  InOEX 

POUTINt  TO  CALCULATE  BFSSEL  MINI  IIO.S  J Sub  ll  F a- 
M Sufl  Ll  OF"  KM  ANO  TmEIk  Of«tVA1|VES 
CALL  BESSLIMMI)  ,KAM) 

J SUM  N Of  KMW  - NOT  PSI  FNC  Of  Ml*  SCATTEPIno  sOLoTIO’* 

SET  IIP  OF  VOM  PI  FNC 

P0*0. 

PlBl. 

P2*3.»A<I> 

VON  TAU  FNC  SET  UP 
TObO. 

tibxiii 

T2b3.*I?.»X(I)*xI1)-I.» 

P bUM  L SOPtR  l SE1  UP 
PLl»SOMT<l.-xn>BA(I>> 

PL2*3.*X ( 1 ) *PLl 

NN  IS  NUMBER  OF  TEPMS  IN  FIELO  * XP ANSI ON 

UO  S L 1 • 1 . NN 

ALbFLOATILI) 

ALL>AL*(AL.|.> 

ALP»AIb*1LI-1)*I2.bAL*1.I/AL/(Al»I.» 

L2*NN»L1 

L3*2*NN«Ll 

L**3*NN«L1 

if  a l.oi.i >ooto  si 


I 


ANi.»  ft  I'm 


J* 


L J 
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TABLE  (I. I).  CONTINUED 


p*p  l 
Tall 
PL*Pl  I 
00 1 0 >>* 

si  If  CLI.uT.llooTO  s<> 

P ap<> 

Tali 
PL»PLi 
00 10  S'* 

NP  P*(i.M./l*l-l.H*i'»l>  *p,?-p  )•(  t.»|./(ai-l.)> 

!•  v 1 | >*(P-p|I-C,».*al-I.I*C  l .-Al  I »*•  « ll  l*P<»*U 
Pl*P<> 

Pi*P 
1 l*Ti 
?i»l 

PL*  li.*l  ./  IAL-1  . > »*»l  1 ) *Pl  «*“PLl«»l  .♦l./IAL-l.M 
PL  I *PL/ 

PL /’•PL 
S'*  CO*?  I Not 

it  NOTH  UOtNUP'  t'KOJCl"*' 

PU  !•*<  1 1 

PLLi»IJ.#»»l>,*«l'-l*>/*. 

PI  ttINCIION  PPO.ItL  TOPS 
PLL  l ■ I • 

PLLi«J.*Attl 

t IPSI  LlOtNIOtt  POLtNOH | AL S 

pll  i «sot»r<i.-*<  i i*«  mi 

PLL<**J.«»tl»»PLLl 
00  S Nla|«.*N 
ANN"?  COAT (Nil 
N*aNN»Nl 
NJ»i*NN*Nl 
NA*J*NN*Nl 

if ini.oT.DoOTo  m 

PLL»PLLI 
0010  bi 

b|  1?  CNl.i.I.ilOOlO  b.l 
PLL*PLLi 
oOlO  bi 

t>  l PLL*  I li.»A«N-I.»**  t I l*PLLi>-l»NN-l  . > ‘PU  l I / ANN 
F IPS?  LtclNOPtS 

h » PLL*  I li.*ANN-l . I*A  ( 1 I *PLL^-»NN*PLL  I ) / UNN-I  . I 

be  PLl  * lANN-l  .»>•»(  I»*PLL*?-PLL  I * U .•!./(  AN  -».n 

PLL l*PLL? 

PLLi*PLL 

Cl  API  ALPnA  SUM  I .N.l  InTTOPAL" 
bi  Cl  1 1>*ALP»ALL#AH(1  I l/MMl ) /PR  (I  >*“1  *PLL**  1 1 » •A«.*«(  i l I 
Cl  t^)a*LP*AHILn*P*PLL#all  I 

Cl  I J I a ALP* I AH 1L 1 * /PR  II) *0Ap IL 1 M • 1 *PLL*a ( 1 1 
Cl («|aALP*AnlLI)*l*PLL*all I 
Cl  lS>aALP*«ANlLH/-M  I»  ♦OAPILU  > • ►’•PLL**  ( 1 I 
Cl  (bl  *A  t CL  II  /KrN*  ALL  *pp  I I I *PL  *P|  L*  I ll/PM  1 >*A-v 
Cl (?)■*  1 1 L 1 ) ApaPLL • A C 1 ? 

Cl  (H)«  (AJILII/PRP'.CAJCLI  I »*l  *PLl  •«  • I I 

Cl  (OlaA.IILl)  *?aPLL*A  111 

Cl  1 10)  ■ CAJCCl  » /KRP*0AJCL1  I l*P»Pl  L»  C I * 

Cl  C I I I *PH  CL  1 I/PR  II  >/PRC  | )*.»R*ALL*P  1 1 • «PL*PLl  •-Hi 
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L).  CONTINUED 


CI<l<*)*->ntLl)  •►'•PLL**  < I > 

Cl  <n)aHH(Lll/<*KI|l  t 

Cl  (U»a-M»L|I«T*PU.»"«1I 
Cl  ( l5l»<**mLl>/PM  I T’OSHIL  I H«p«PLL#a<  i I 
CALL  AOO(H.L?.L3«LA.N1 

C ArflO  »&.9JILl«NI.ALP«Aj.AN.rtH,.lA|.  .4rt.U8M.X  <1  > .P.  I .PLL.PL 

9>  FORMAT  <SX<<TI10<2  15>  t£E20.MI  ,/J  <^x,  st  io.Ml  ./,?<“>«  ./>-  n . - I . 
ISX.SL20.8) 

5 COM  1NOF 
A COMTINOC 
C CALL  API 

CALL  CGAU55INNN.1  .•..NNN.U.  l«.t.  >» 

«P I IF  (6.V0H- 
90  FO^**Ar</5A»»<!0.81 
CALL  AP3 
9h  CO'iT|NUC 
99  STOP 
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TABLE  (1.2).  LISTING  OF  ROUTINE  START 


I" 


Is 


in 


2s 


IN 


SUsROUTtNt  SlAwT 

C0MFFON/0NE/XT20)  .*(20*  .AM.ALAM.SM.  IN.  4.  < JS) 

coHM.)N/TH^tt/A««<;oi  .ajooi  .««(^'ii . Ajum  .iuh (->.!> . ihi  >n 

COMPLEX  AJ.Art»Mrt.o*j.L)AM.I.Mfi 

COMMON/F  OUH/NT  .NP  . KM  I N • l)R  .SU.SP.AXX  ( ls>  .nxx 
COMPLEX  AM 

NAMtLlsT/lNlO/AM.AL  AM.WMIN.WMAX.SH.3H.I-N.NT  ,n  M.(W 

c su«houtinl  start 

C SUiNOUTINt  It)  SET  '.M  NfCESSARY  ->A~  I«1  lr«s  F Ox  *A  I N 

C %H»**«4###**4v***«***  .«••««  ••«**« 

C X IS  ARRAY  of  GAUSN-LEGtNOct  COMF1*  I HI  I • VALU'N 

C ■ IS  ARRAY  OF  CORN*  sPONOlNi. 

C 20  POINT  G-L  OllAORAloPt 

C Exi-ECTtO  URPtH  Lt*M  I TMRRtFORt  1 s 

C AM  Is  COMF'LFA  !NOF>  OF  RFF-aCTT  'n  RarTICLF 

c an  is  array  of  surface  coffficifnt- 

C HUM  rtt  IN  FORM  OF  ttROTM  LtGtNoPF  POLYS 

C ALAM  Is  rAvELENG Tm  IN  ONI  I S OF  M|f  'ONN 

C NM  Is  NllMbt  R OF  Tt-MS  In  SuRFACF  F«RANs|ON 

C NN  Is  NUMOt  R OF  TF  RMS  TO  MF.  C0NN|O*  RfcU  IN  raPAlSlO 

C 

UATA(XI!>.1*1.10>/.0  76N?6S/II../?7  7<«SmM  l . . 3 7 J fami-h  7 , .M  0 -s  70 n I R . 
1 .6T6uSJN«0  7,.T4ti33l'*0tu...F»  (M  1 -»R7  I - . .H  I 22T«.*202<  .s  JR  7 \ *.2  7/. 
?.RH312bS,»M/ 

OAIA(*(t»*l*|.10»/.ls27S33"7I..|«.R>  72H**,...  . I «•/(  RM  v J.  . I 31  • 

1 .1  K»l4*S3lN,.  10 IV  Till  ISrt.  .OR  3?7n7aI  -.  .OnaM/aam  *.  I4/9-. 

2.01 7M<*00  71/ 

00  1 1*11*20 
11*21-1 
X<l)*-A<lll 
1 « < t > *M  < 1 1 > 

REAOIS.INlOl 
■HIT. (b.INTol 

OR=<RMAX-H'MN> /FLOAT  (NT  -1  > 

HE  TOON 
ENi' 


■■I  I 
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TABLE  (1.3).  LISTING  OF  ROUTINE  AP2 


SU'MIUlINt  AP<» 

CUMMuM/KNt/H  tiin  I »*  1^0)  .AM.aCAM.NM.nN.AM  15) 

Common/  ( 40/  M HO  «tt  1 ) .Cl  (IS)  •*! 

COMMON/)  OOK/NT.NP.-MIN.iIP.SU.SP.AXXI  IS)  .NXX 
DJUnSION  mm  uni  I . P/T  ( 101)  .PHI |0  1 1 .PUT  (101 > 

COMPlF/  A.Al.CI 
CO-'PlEA  AM.  ON 
CO*PltA  SO'.k.SlAO 
If  INp.of  .1  M.olo  SP 
ST  =0  . 
t T*0, 

ft*o. 

ou  s » t«i . lot 

PI i ( i i *o. 

P3 ' t|)*n. 

PM  r ( 1 I *11  . 

Si  P?I  I 1 1*11. 

s ? wpi  n it.,')/) 

fP  = F »p<-(  ( «U0&  I AN  < 1 ) )-ALOO|SPI  > /<A|  00  ISO)  )**2/P.)  1/ANI  I > 
F 1*F T.Fo 
NNo*«  *NN 
1)0  i 1*|,,.NN 
*P|  Tt  (*i,PU»  1 • A ( 1 • 1 ) 
o<  io-Miir<sx.iio.pE^u.a) 
p C0>  iT  | NOF 
0/  F O-MAT ( 1 H 1 ) 

1)0  3 1 * | . NN 
N*P*NN»| 

NA* }*NN«I 
fN*f  COAT  <l> 

liN*(-Al>**t|M>*fN»CFN.l.>/(P.*FN«|.l 
. A ( I.  1 ),*A  (N.  I ) *(.N 
l A(SA.l)*A(MA.l)»OS 
■PI  11  lo.VP) 
nnnp*/*nn* I 

00  N i*NNNP,NNN 
■PIT) (m.PU) I.A(l.l) 

CONTINUE 

NNNX*0 

s*u. 

1*11, 

00  S I*|  .Nl, 

**  ■ ♦ I • • 

I.A*  3»NN.  I 
AM*F|  OAT ( I ) 

A A*«?  . • Am  . 1 . 

f *H  APS  (A(N.l)  )**P.CAHS<  AINA.  !))••/) /Am 

NNFJXa  | 

S*--.F  *am*A« 

t«-  • IPF.AL  I AIN.  1 ) ) «>-tAL  ( AIN*,  1 ) ) )*A« 

If  If  .Lt.l.t  -1AI0OII)  SI 
m CO  j I 1 NOF 

Ml  F AI.*ALAM*AI  AM//./,).  lAlSPi? 

S»''*F  A(. 
t*>  *f AC 
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SHIS  PAGE  IS  BEST  QUALITY  PRACTICABUI 
OOP!  VUtiJUlSlUM  w lADC  


TABLE  (1.3).  CONTINUED 
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TABLE  (1.4). 


LISTING  OF  ROUTINE  ADO 
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TABLE  (1.5).  LISTING  OF  ROUTINE  CGAUSS 
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TABLE  (1.6).  LISTING  OF  ROUTINE  BESSL 
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APPENDIX  II 


SPHERICAL  BESSEL  FUNCTION  GENERATION;  NUMERICAL  ASPECTS 


One  of  the  difficulties  experienced  in  this  investigation  was 


the  calculation  of  the  radial  functions  and  P»(x),  Eqs.  (2.45) 


and  (2.47),  and  their  derivatives.  This  difficulty  is  not  new,  and 


in  the  early  days  of  the  implementation  of  the  Mie  solution,  it  was 


their  calculation  (Deirmendjian,  1964).  Since  that  time,  however 


advances  in  the  construction  of  compilers  has  permitted  more  fundamental 


calculations  of  the  radial  functions  if  sufficient  safeguards  are  taken 


Because  of  the  special  dependency  of  this  investigation  on  calculating 


the  values  of  these  radial  functions,  this  appendix  describes  the 


numerical  considerations  of  calculating  these  quantities 


The  radial  function  i|>»(x)  and  P»(x)  are  related  to  the  spherical 


bessel  functions  and  spherical  hankel  functions  of  the  first  kind  by 


where  V »(x)  » spherical  Neuman  function.  The  first  two  tp»(x)  and 


Pq(x)  ■ sin(x)  + -tcos(x) 


(II. 5) 


Pj<„>  - - CO.(x)  + + >1"W) 

Additionally,  the  ^(x),  j/^(x),  and  h^(x)  satisfy  the  same 


(II. 6) 


recursion  relation 


f»(x) 


f£+1(x)  - (2t+l)~  - ft_lM 


(II. 7) 


and  derivative  formula 


- f£-i(x) 


(II. 8) 


where  f^(x)  may  be  any  of  the  three  functions.  Additionally,  the 
j^(x)  and  ij^(x)  have  a cross  product  relation 


(*)(/£_!  (x)  - y^i(x)^(x)  - — j 


(II. 9) 


The  recursion  relation  for  (x ) or  P^(x)  may  be  obtained 


from  Eq.  (II. 7)  by  multiplying  by  x to  yield, 

♦/(x) 

*£+1(x)  * (2£+l)- (x) 


and  the  derivative  relation  may  be  similarly  obtained  as 


(II. 10) 


dt  ~ *t-iM  ' x Vx) 


(II. ID 


In  principle,  Eqs.  (II. 10),  (II. 11)  and  (II. 3)  - (II. 6)  may  be 
used  to  calculate  the  value  of  ^(x)  and  P^(x)  for  any  x and  ( by 
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forward  recursion.  In  practice,  however,  the  forward  recursion  tends 
to  diverge  numerically  when  ]x|<l,  causing  problems  with  overflow  as 
well  as  inaccuracy  due  to  roundoff  errors. 

To  circumvent  this  difficulty,  the  calculational  procedure  of 
Miller  may  be  used  (Miller,  1950).  This  procedure  operates  as  a back- 
wards recursion.  The  algorithm  operates  in  the  following  manner: 
some  maximum  value  L for  l Is  established.  An  upper  limit  for  the 
algorithm  L*  at  least  twice  L in  value  is  established.  The  larger  L* 
is,  the  greater  the  accuracy  of  the  method.  A series  of  factors 
are  calculated  using  Eq.  (II. 10)  in  a backwards  recursion  subject  to 


the  conditions 


F - 1 
L*  1 


(11.12) 


F * . 0 
L*+l  U 


(11.13) 


The  may  then  be  calculated  from  the  F^  by 

**(«)  - Vx>VF0 


(11.14) 


where  iJ>o(x)  *s  defined  by  Eq.  (II. 3) 

It  may  be  noted  that  this  method  also  suffers  from  divergence 
problems  in  a like  manner  as  the  forward  recursion  technique.  Fortu- 
nately, this  does  not  occur  within  the  bounds  of  the  Mie  solution  or 
this  investigation  for  |x|<  18.  In  practice,  the  two  methods,  forward 
and  backward  recursion,  overlap  so  that  values  of  i^(x)  may  be  calcu- 
lated for  any  except  very  small  values  of  x where  Rayleigh  or  Rayleigh- 
Gans  theory  is  valid. 


fflk 
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I 

The  P^(x)  could  also  be  calculated  using  Miller's  method,  but 
are  not  because  It  is  easier  to  rewrite  Eqs.  (II. 2)  and  (II. 9)  as 

P^(x)  - ^(x)  - ■cn^(x)  (11.15)  | 

and 

X»^(x)n^  j(x)-l 

Vx)  " <n-16> 

where:  n^(x)  - X^x),  and 
nQ(x)  - -cos(x). 

The  n^(x)  may  then  be  calculated  from  the4^(x),  and  the  p^(x)  from 
the  ^(x)  and  the  n^(x). 

The  routine  BESSEL  uses  this  recipe  to  calculate  the  values  of 
the  radial  functions: 


1.) 

Calculate 

the  F^,. 

»•  V 

2.) 

Calculate 

the  i^(x) 

from  the  F^. 

3.) 

Calculate 

the  P ^(x) 

(via  the  n^(x))  from  the  if^(x). 

4.) 

Calculate 

V 

the  \|>^(x) 

and  the  p ^ (x)  from  the  t^»(x)  and 

theP^x)  via  Eq.  (II. 11). 

Based  on  a survey  of  the  literature.  It  appears  that  Miller's 
method  had  previously  been  validated  only  for  real  value  of  x.  During 
this  investigation,  therefore,  its  validity  for  complex  values  of  x 
was  demonstrated  through  the  numerical  validation  of  the  entire  code. 


APPENDIX  III 

REDUCTION  OF  THE  NONSPHERICAL  FORMALISM  TO  MIE  THEORY 
FOR  THE  SPHERICAL  PARTICLE  CASE 
As  a matter  of  course,  the  nonspherical  formalism  developed 
in  Chapter  VI  must  reduce  to  Mie  theory  when  the  particle  is 
spherical.  This  appendix  will  demonstrate  that  reduction. 

In  terms  of  Eq.  (6.1),  the  radius  of  a spherical  particle  reduces  to 


a(9)  = S P (cos0). 
o o 


(III.l) 


and  has  angular  derivative 

0. 


3a(6) 

36 


(III. 2) 


Additionally,  the  radial  functions  ^ and  and  their  derivatives 
are  independent  of  0 , so  that  Eqs.  (6.43)  - (6.57)  reduce  to 


ku\  i * 0 


^(ka) 

V’,2  * yi  \tr 


4^!  (ka) 


'U\  3 - H ka 


^(ka) 


'U\U  ~ H kT 


hit' 


l2W 


(III. 3) 


(III. 4) 


(III. 5) 


(III. 6) 
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ic 

- 


p^1  (ka) 


AC£’,15  " ka  *2 U' 


where 


U 1 


?t  ?r 


sin(0) 


dcosO 


f3Pt’  1 

l2 it'  m J 36  PV  dcose 


Equation  (6.63)  may  then  be  written  as 


i^,,(kma)^  \iy(kma)  ty^i(ka) 

kma  iti  * l~»a. 


kma  I2lt'  ~ka  1 lit' 


- V (kma)T  -^<kma)  tp/(ka) 

kma  2(Ji ' kma  ^ ^ p p i 


,(kma)  _t4»<.(kma) 

kma  ^2ii’  kma  ltt'  ka 


P^ (ka) 


I T p p I 


(III. 17) 


(III. 18) 


(III. 19) 


-tp^1(ka) 


I^ppi 


p^,1  (ka) 


tp^1 (ka) 


ipi^(kma)  (J^(kma)  (ka)  p^l'(ka) 

*iU'  ]2it'  ki  * ki  \u' 


(III. 20) 


ana  Eq.  (6.64)  as 


' 1 
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ittilttAiBIlitoiriiiiin'rriilr'i 


(III. 21) 


/ \ 
) Y^> ( (ka)  1 ^ 'J,£*  (^a)  1 1££ 1 I 

I f*J^  (ka)  ( 2f £ 1 (ka)  1 2 • 1 /^® 

—j^Y^ (ka)  I *Pg  (ka)  1 » 1 /ka 
*Ty^  [ (ka)  I (ka)  ] /ka 


To  simplify  the  demonstration,  two  restrictions  may  be  made 
without  compromise:  first,  that  the  particle  is  small  enough  that 
only  the  first  term  ( - 1 need  be  retained,  and  second,  that  V he 
chosen  such  that 


1itr  ’ 1 


(111.22) 


l2W  ’ 0 


This  reduces  Eq.  (111.20)  to 


(III. 23) 


l-  * 


(kma)  0 

(kma) 

— 0 

m 

0 -ip  *(ka) 


and  Eq.  (III. ill)  to 


(III. 25) 


It  may  be  noted  that  Eq.  (III. 24)  is  now  separated  In  terms  of 
pairing  a^  with  c^t  and  b^  with  d.^.  If  Eq,  (III. 24)  and  (111.25)  are 
combined,  the  result  is 


I 1 • 1 

(kma)  - c.p.  (ka)  - y (ka) 


(III. 2b) 


Equation  (III. 26)  may  be  rewritten  as 


(ka)]/4»'  (kma) 


and  combined  with  Eq.  (III. 29)  to  yield 


(ka)<l>  (kma)  - (kma)i^ 


*%*r 


i 


I 

i 

: 

\ 

i 


I 

{ 

i 

I 

I 


4-CjJpj  (ka)^  (kma)  - mp*  (ka)<^  (kma) ] - 


tYj  [mi^  (kma)^ ' (ka)  - (kma)<l/(ka)] 


[mij^(kma)t|^  '(ka)  - '(kma)^  (ka)  ] 

c Yj  


[p11(ka)^1 '(kma)  - mp^  (ka)^  (kma)  ] 


(III. 33) 


(III. 34) 


i 


which  Is  Identical  to  Eq.  (3.17)  for  l * 1. 

In  a like  manner  Eqs.  (III. 27)  and  (III. 28)  may  be  combined  to 
yield  Eq.  (3.18)  for  ( - 1. 


- . 1 


r 


„ 


sr a as 


&,  ' 
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APPENDIX  IV 

TWO-DIMENSIONAL  RADIATIVE  TRANSFER  CODE 

The  computer  code  described  in  this  appendix  implements  the 
algorithm  of  Chapter  X.  Calculations  performed  using  this  code  are 
presented  in  Chapter  XI.  The  version  of  the  code  presented  here  is 
for  a constant  extinction  coefficient. 

The  computer  memory  requirement  for  this  code  is  approximately 
1.25  x 8^  words.  This  requirement  is  not  decreased  appreciably  by 
compilation  using  the  SCOPE  3.4.2  compiler.  Run  time  for  a nominal 
run  in  the  near-ir  requiring  24  orders  of  scattering  to  satisfy  the 
1%  self-consistency  requirement  is  about  five  seconds.  The  same 
calculation  for  varying  medium  would  require  about  twice  as  much 
time. 

The  code  consists  of  a driver  routine  and  five  subroutines. 
These  routines  are  described  individually  below,  and  are  listed  in 
Tables  (IV-1)  - (IV. 6). 

RT2  (Table  (IV. 1)) 


This  is  the  driver  routine.  It  inputs  the  extinction  coeffi- 
cient, the  albedo  of  single  scattering,  and  the  phase  function  for 
directions  i = 1 to  5,  sets  the  scale  of  the  calculations,  and 
normalizes  the  phase  function.  The  calculation  is  performed  for 
incremental  ranges  or  optical  depths.  Routine  LSTART  is  called  to 
generate  the  zeroth  order  of  scattering  source  functions.  Routines 
UPDOWN,  ACROSS,  and  ACCOUNT  are  called  for  each  order  of  scattering 
to  calculate  the  source  functions.  Finally,  CTRST  is  called  to 
calculate  to  modulation  contrast  for  the  total  Intensities. 
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LSTART  (Table  (IV. 2)) 


This  subroutine  uses  the  Incident  Illumination  beam  profile, 
the  extinction  coefficient,  and  the  phase  function  to  calculate  the 
zeroth  order  of  scattering  source  functions.  This  routine  is  specific 
for  beam  type  sources  with  light  propagating  only  along  the  4.  * 1 
direction.  More  general  incident  illuminations  are  not  considered  by 
this  routine  because  only  this  type  of  incident  illumination  was 
needed  in  this  investigation. 

UPDOWN  (Table  (IV. 3)) 

Subroutine  UPDOWN  calculates  the  nth  order  of  scattering  source 
function  contributions  for  directions  4,  * 1,  3,  5,  7 from  the  n-lth 
order  of  scattering  source  functions.  Additionally,  it  calculates  the 
nth  order  of  scattering  intensities  for  directions  4.  * 1 and  5 at  the 
edges  of  the  medium.  These  are  the  forward  scattered  and  back- 
scattered  intensities,  respectively. 

ACROSS  (Table  (IV. 4)) 

This  routine  performs  the  same  source  function  calculations  as 
UPDOWN  except  for  directions  4.  m 2,  4,  6,  8. 

ACCOUNT  (Table  (IV. 5)) 

Surboutine  ACCOUNT  transfers  the  nth  order  of  scattering  source 
function  array  into  the  n-lth  array  in  preparation  for  the  next  order 
of  scattering  calculation.  Additionally,  it  zeroes  out  the  nth  order 
array. 
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Integrals  used  in  non-spherical  particle  scattering 
solution. 

Field  expansion  coefficient  in  Rayleigh  Hypothesis. 

Reilly  integral. 

Arbitrary  vector  function. 

Matrix  in  least-squares  technique. 

Transpose  of  matrix. A. 

Inverse  of  matrix. A+A. 
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nth  order  of  scattering  source  function. 

J^(x) 

nth  order  of  scattering  source  function  at  ^ along  J^. 

Spherical  Bessel  function  of  the  1st  kind  of  order  * . 

I to' 

Propagation  direction  unit  vectors. 

h'V 

ith  propagation  direction  unit  vectors. 

k 

Wave  number. 
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H 

■L - 

• * [ 

1,1' 

t 

In 

m 

ml’mo 

N(^,a,t) 

n 

n 

n(r) ,n(a) 
P(^,t,0),P(e) 
P^(cos0) 
P^(cos0) 

Pj  Cfc.t*0) 

V 

Q(^.t,0) 

R(r) 

*»  t* 

I’l 

fyjn 


Summation  indices. 

Path  length. 

Natural  logarithm. 

Refractive  index,  summation  index. 

Refractive  indices  inside  and  outside  particle. 
Particle  distribution  function. 

Surface  normal  vector. 

Unit  surface  normal  vector. 

Index  of  order  of  scattering,  summation  index. 

Particle  size  distribution  function. 

Unpolarized  phase  function. 

Legendre  polynomial  of  order  t. 

mth  associated  Legendre  polynomial  of  order  l. 

Scattering  phase  function  corresponding  to  jth  Stokes 
vector  component. 

Polarized  phase  function  used  in  radiative  transfer. 

Weighted  quadrature  approximation  phase  function  value 

for  scattering  from  V.  to  k .. 

i J 

Modified  unpolarized  phase  function. 

Reilly's  unknown  function. 

Radial  unit  vector. 

Position  vector. 

Position  of  the  l,m  point  in  L x M array. 
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s,s' 

S1(0),S2(0) 

"s" 

Sn 

sin 


U 

s 

u 

V 

«y» 

x‘Xo 

x(0) 

y^(x) 


ai 

abn 

e<*.o 

Y 


Particle  radius,  radial  coordinate. 

General  solution  of  the  Helmholtz  equation. 

Angular  scattering  function. 

Superscript  indicating  scattered  wave  or  field. 

Expansion  coefficient  in  a(0,4>). 

Sine  function. 

Time  variable. 

Stokes  vector  component. 

Scalar  scattered  field.  * v • 

Stokes  vector  component. 

Superscript  indicating  interior  ("within")  wave  or  field. 
Mie  size  parameter  («2Tta/X). 

t 

Irregular  particle  shape  function. 

Spherical  Bessel  function  of  the  2nc*  kind  of  order  t. 

Unit  vector  in  the  z direction,  identical  with  e^. 
Cartesian  coordinate. 

Angular  scattering  function  expansion  coefficient. 

Extinction  coefficient  ar 

Extinction  coefficient  at  £ at  time  t. 

Angular  scattering  function  expansion  coefficient. 

Scattering  coefficient  at  £ at  time  t 

Parameter  for  Deirmendjian  particle  size  distribution 
function. 


6 


I 7 M 


Incident  Debye  potential  expansion  coefficient. 

Parameter  for  Deirmendjian  particle  size  distribution 
function,  phase  difference. 

Dielectric  function  (constant)  total  fitting  error. 


PI 


a 


I 

j 

eil,€U 

Fitting  errors  in  least-squares  technique. 

n 

Exponent  coefficient  in  Junge  particle  size  distribution 
function. 

n^(x) 

Radial  function,  * x^(x). 

n^(0) 

Angular  function  used  in  S^(0)  and  82(6). 

K 

Real  part  of  m. 

1 

X 

Wavelength . 

u 

Permeability. 

« 

u 

Cos(6),  • z. 

V 

Negative  imaginary  part  of  m. 

j 

U) 

Angular  frequency. 

• 

UJ* 

Effective  albedo. 

1 

w(j£,t),U) 

Albedo . 

! 

• 

4> 

Magnetic  scalar  potential. 

i 

Scalar  Helmholtz  equation  solution  of  kellly. 

■ 

4>  (r,0,$) 
n 

Symbolic  linearly  independent  function  cf  least-squares 
technique . 

<t> 

Azimuthal  angle. 

1 

l 

Electric  Hertz  vector. 

1 
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£ 

* 

i 

h 

U 


\(r,e,4.) 

^(x) 

p 

1 2 
p£(x),  P^(x) 

$ 


u 


0 (a,A,m) 

ex 

oi(a,A,m,0) 
o (a,A,m) 

SC 

Osc(a,A,m,0) 

V0) 

0 

0 

<r> 


CD 


Electric  Debye  potential,  ratio  of  circle  circumference 
to  diameter. 

Electric  scalar  potential. 

Symbolic  linearly  independent  function  of  least- 
squares  technique. 

Radial  function,  * Xj^(x). 

Charge  density. 

1 2 

Radial  function,  ■ xh^(x)  and  xfi^(x) 

Magnetic  Hertz  vector. 

Standard  deviation  used  in  log-normal  particle  size 
distribution  function,  conductivity,  magnetic  Debye 
potential. 

Extinction  cross  section. 

Stokes  vector  component  i differential  cross  section. 
Total  scattering  cross  section. 

Unpolarized  differential  scattering  cross  section. 
Angular  function  used  in  S^(6)  and  82(6). 

Polar  angle  and  scattering  angle. 

Unit  polar  vector. 

Gradient  operator. 

Boundary  condition  E function. 

Boundary  condition  function. 

Average  particle  size  used  in  log-normal  particle  size 
distribution  function. 

Infinity. 
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